Use of mathematical modeling to study pressure regimes in normal and Fontan blood flow circulations

We develop two mathematical lumped parameter models for blood pressure distribution in the Fontan blood flow circulation: an ODE based spatially homogeneous model and a PDE based spatially inhomogeneous model. We present numerical simulations of the cardiac pressure-volume cycle and study the effect of pulmonary resistance on cardiac output. We analyze solutions of two initial-boundary value problems for a non-linear parabolic partial differential equation (PDE models) with switching in the time dynamic boundary conditions which model blood pressure distribution in the cardiovascular system with and without Fontan surgery. We also obtain necessary conditions for parameter values of the PDE models for existence and uniqueness of non-negative bounded periodic solutions.


Introduction to Fontan surgical procedure
In a normal biventricular heart, the systemic and pulmonary blood circulations are in series and each one is supported by a ventricle. The Fontan surgical procedure is applied to a malformed heart that is characterized by the presence of only one (left or right) functional ventricular chamber, and for which a biventricular repair is not possible [MMM + 13]. The Fontan procedure was first introduced in 1968 and involved routing systemic venous blood flow to the pulmonary arteries, without use of an additional ventricular pumping chamber. In its modern form, the Fontan procedure involves creating an anastomosis (surgical connection) between the superior vena cava (the vein responsible for carrying blood from the upper body to the heart) and the left and right pulmonary arteries, and between the inferior vena cava (the vein responsible for carrying blood from the lower body to the heart) and the left and right pulmonary arteries, with the connection between the inferior vena cava and the pulmonary arteries being made using a synthetic graft called an extracardiac conduit. The resulting connection between these four vessels is referred to as the total cavopulmonary connection. The Fontan circulation creates an unusual state in which the force driving the pulmonary blood flow is solely a residue (in the systemic venous pressure) of the only functional ventricular chamber's contractile force. This single ventricle propels blood flow through the systemic arteries and capillaries, with the systemic venous return passively entering the pulmonary circulation. The ventricle is doing nearly twice the expected amount of work because it has to pump blood to both the body and the lungs.
Fontan surgery is an extraordinary story of success in that it has allowed a generation of newborn babies with the most severe forms of congenital heart disease to survive into adulthood (estimated prevalence of approximately 1 per 3000 births) [KPM07]. Though life-saving, a univentricular Fontan circulation does not, however, reproduce biventricular physiology and has been considered abnormal in the sense that systemic venous hypertension (mean pressure > 10 mm Hg) occurs simultaneously with pulmonary arterial hypotension (mean pressure < 15 mm Hg) [DL05]. In patients with Fontan physiology, life expectancy remains far below projected age-and sex-matched normative values. Patients with Fontan procedures most frequently die from heart failure or from thromboemboli [KFM + 08]. The incidence of thromboembolic deaths rises sharply 15 years after Fontan surgery [VBC + 09, dIC + 07]. At the same time heart failure deaths are very uncommon during the first 10 years after Fontan surgery, with a steady decline in survival thereafter. Associated factors include protein-losing enteropathy, single right ventricle morphology, and increased Fontan pressures [MMM + 13].
Late Fontan failure might progress gradually over years with an absence of overt symptoms. Fontan patients have lived with less than ideal cardiac output their entire lives and might not recognize decline in functional status until deterioration is significantly advanced. In the medical literature, failure of the Fontan circulation is divided into 3 main categories: ventricular dysfunction, systemic complications of Fontan physiology, and chronic Fontan failure [GSRR11]. In a cross-sectional analysis of 546 children with Fontan procedures, 27% had abnormal ventricular ejection fractions and 72% had diastolic dysfunction. The prevalence of systolic and diastolic ventricular dysfunction continues to increase in adulthood [PVS + 02, EFG + 03]. Over the past 5 − 10 years, a number of studies have described the effect of Fontan physiology on the liver. Hepatic venous pressure after Fontan surgery might be 3 − 4 times higher than normal, with levels commensurate with congestive heart failure in adults. Implications of this elevated hepatic venous pressure over the long term remain to be fully understood [KSH + 08, CMS + 08]. Fontan physiology is characterized by progressively decreasing cardiac output and increasing central venous pressure over time. The average peak oxygen consumption ranges from 19 to 28 mL/kg per minute, or 50% − 60% of predicted values [FMK + 10, PMC + 08]. In the third decade of life, hospitalization rates and symptoms increase significantly [DGD + 10]. In the future, a "Fontan pump" or a "Fontan assist device" might usher in a new era of ventricular replacement. Although investigations have begun, such devices still appear to remain many years away [RCF + 10, BKCT09]. Recently, the use of the right ventricular Impella has been described to direct flow into the pulmonary artery, resulting in a modest reduction in central venous pressure [HFTM + 12].

Computational fluid dynamics
Computational fluid dynamics is a powerful tool that can be used to gain insight into the local blood flow dynamics in the Fontan circulation. These simulations are used to model the detailed 3-D hemodynamics of a particular region in the body, such as the total cavopulmonary connection, rather than the complete circulation. A simplified three-dimensional model was used [THZ96] to simulate the local fluid dynamics for different designs of the total cavopulmonary connection, allowing a quantitative evaluation of the dissipated energy in each of the examined configurations. The authors show that, from a comparative point of view, the energetic losses can be greatly reduced if a proper hydraulic design of the connection is adopted, which also allows control of the blood flow distribution into the lungs. Under the assumptions that vessel walls are completely rigid (according to surgical reports, the vessel diameter change per cardiac cycle is around 5 − 10% in most of the major arteries) and all vessels are symmetric, numerical simulations of blood flow to the lung after a surgical Fontan procedure are described in [DDLP + 96]. The authors simulated the full nonlinear Navier-Stokes equations using a streamline finite element method. They analyzed the blood flow dynamics for different values for the offset between the superior vena cava and the inferior vena cava anastomoses and concluded that the optimal distance for the offset is about 7 mm. It was shown in [BHB + 09] that wall flexibility can play an important role in determining quantities of hemodynamic interest in the Fontan connection. However, [TWF + 17] recently showed that while fluid-structure interaction effects are important for instantaneous quantities of interest in the Fontan circulation, they have a negligible impact on time-averaged values.
According to the article [DeG08], the main quantities of importance in modeling the Fontan procedure are: • Vessel diameters and flow rates representative of the range seen in the patient group under study including resting and exercise states • Vessel sizes and flow rates matched appropriately • Compliant vessels, accurate modeling of surgical anastomosis sites, and surgical material used (unless proven unnecessary) • Unsteady flow • Effects of respiration • Correctly shaped vessel anatomy Two different types of boundary conditions, time-averaged and pulsatile, were analyzed in [WTT + 16]. The authors derive a patient-specific sensitivity criterion which provides a guideline for determining when time-averaged boundary conditions can be used to save computational time.
Recent advances in imaging methods and patient-specific modeling now reveal increasingly detailed information about blood flow patterns in healthy and diseased Fontan states. Building on these tools, there is now an opportunity to couple blood flow simulation with optimization algorithms to improve the design of surgeries and devices, incorporating more information about the flow physics in the design process to augment current medical knowledge. To do so, there is a need for efficient optimization tools that are appropriate for unsteady fluid mechanics problems, particularly for the optimization of complex patientspecific models in the presence of uncertainty. The state of the art in optimization tools for virtual surgery, device design, and model parameter identification in cardiovascular flow and mechanobiology applications are reviewed in [YFM10]. In this work, the authors perform optimization on a model Y-graft design problem. This work represents the first use of formal design optimization methods for the Fontan surgery, and also demonstrates the applicability of the optimization framework on a pulsatile flow problem with multiple design parameters and constraints.

Lumped parameter models
While computational fluid dynamics models can be used to calculate detailed three-dimensional blood flow in the total cavopulmonary connection, the computational costs of this approach prevent it from being used to simulate the entire circulatory system. Because Fontan failure is a systemic problem, reduced order methods, such as lumped parameter models, can be used to study this phenomenon. Lumped parameter models are based on the analogy between fluid flow and electric circuits. In addition to modelling the entire circulatory system, these models can also be used to generate appropriate upstream and downstream boundary conditions for computational fluid dynamics simulations.
A lumped parameter model of the Fontan circulation was used by [TCT + 11] to generate boundary conditions for a computational fluid dynamics model used to design a Fontan assist device. In a study by [LSK + 14], lumped parameter models of the Fontan circulation and the normal circulation were compared to determine differences between the two circulations in the regulations of cardiac output and central venous pressures. In a study by [KPM + 14], a lumped parameter model was used to study the Fontan circulation under exercise conditions.
The objective of the present study is to develop lumped parameter models of the Fontan circulation with the goal of understanding the systematic changes that occur during Fontan failure. The outline of this paper is as follows. Section 3 describes an ODE model of the Fontan circulation and presents some basic results from this model, showing its consistency with physiological behaviour. Section 4 describes PDE modeling of the normal and Fontan circulations, which is an extension of the ODE approach that allows for spatially inhomogeneous variation of the model parameters. Section 5 determines range values of parameters for PDE models from the Section 3 such that a unique non-nenegative solution exists. We also show in the appendix how to construct super-and sub-solutions for PDE models. 3 Spatially homogeneous ODE model of blood pressure distribution for the Fontan circulation A simple model of the Fontan circulation can be based on an electric circuit approach. This model consists of five compartments: the heart, the arterial system, the capillary system, the venous system, and the pulmonary system (lungs). For the Fontan circulation, all compartments are connected in series around a single loop. In particular the flow must pass sequentially through the arteries, capillaries, veins and lungs, before returning to the heart, as shown in Fig. 1. We model the capillary and pulmonary systems as linear resistance vessels. That is, we assume that the pressure drop across the vessel is proportional to the flow through the vessel, with a constant of proportionality called the resistance, labeled R c and R p , for the capillary and pulmonary systems, respectively. We assume these vessels have no compliance. We model the arterial and venous systems as linear compliance vessels, in which the volume of the vessel is proportional to the pressure in the vessel, with constant of proportionality called the compliance. We allow for compliance vessels to have resistance, which is modeled in the same way as for the resistance vessels. The compliances of the arterial and venous systems are labeled C a and C v , respectively, while their resistances are R a and R v . The heart is considered to be a linear compliance vessel with different compliances depending on whether it is relaxed (in diastole) or contracted (in systole). As shown in Fig. 1, the heart has compliance C d in diastole (0 ≤ t < 0.7) and compliance C s in systole (0.7 ≤ t < 1). We assume that all vessels contain some basal volume at zero pressure, which are denoted as V 0 i . The variables in the system are the volumes V a , V v and V h of the compliance vessels (the arterial system, the venous system and the heart), and the pressures P A , P a , P v and P pv at different points along the loop; see Fig. 1. The parameters of the system are the resistances, compliances and basal volumes of all the vessels as well as the total blood volume V T . Estimates for all the parameters can be found by measurement on individuals. For the results that follow, the values of the parameters that we use have been taken from the literature and are given in Table 1. Parameter values: V in units of L, C in units of L/mm Hg and R in units of mm Hg· min/L. It should be noted that V s was chosen to be negative to achieve a reasonable value of systolic heart compliance C s . In practice, the systolic pressure is never a small value or 0, so the volume in the heart never reaches the negative value V s = −0.5 L.
We consider the compliance of the heart to be a piecewise constant function, with value C s in systole and value C d in diastole. To ensure appropriate directionality of the forcing, we assume that there are "perfect" valves where the pulmonary vein enters the heart and where the aorta leaves the heart. Anatomically, in Fontan circulation, the pulmonary vein enters a common atrium which is separated from the single ventricle by an atrioventricular valve. Depending on whether the patient has a functioning right or left ventricle, this valve is either the tricuspid valve or the mitral valve. Similarly, the single ventricle is separated from the aorta by a semilunar valve, which is either the pulmonary valve or the aortic valve. For the purposes of these models, we are neglecting the common atrium by considering it an extension of the pulmonary vein, and we refer to the valve separating the pulmonary vein from the heart as the "pulmonary vein valve" and the valve separating the heart from the aorta as the "aortic valve". In our model, during diastole, the pulmonary vein valve is open and the aortic valve is closed, allowing flow into the heart, and during systole, the pulmonary vein valve is closed while the aortic valve is open allowing flow out of the heart and into the arteries. We assume that the valves are perfect in the sense that they open and close instantaneously and in synchrony, and that they do not allow back flow, regardless of pressure differences or flow characteristics. These assumptions lead to a cardiac pressure-volume cycle as presented in Fig. 2. We use conservation of volume in each compartment to set up the dynamic equations. In particular, we must have the rate of change in volume of a compartment equal to the difference of the flow in and flow out. Due to the assumptions on the valves, there are discontinuities in the variables and their derivatives as the valves switch. Thus, it is convenient to distinguish the time intervals in which the heart is in systole and in which it is in diastole. As such, there is one set of model equations that is valid in diastole, while another set of equations is valid in systole.
The three flow rates are the arterial, capillary and pulmonary flow rates, defined as The three compliance volumes are the arterial, venous and heart volumes, defined as Consequently, the dynamic equations are, during inflow (diastole), so that Q a = 0 and P A = P a , and during outflow (systole) so that Q p = 0 and P v = P pv . In addition, total volume must always be conserved, so that It is possible to simplify these equations by using the conservation law (3.6). In particular, during inflow (diastole), and during outflow (systole) (3.8) In addition to diastole and systole, the cardiac cycle consists of two isovolumetric phases, during which both heart valves are closed, and the heart undergoes a change in pressure in response to a change in its shape, while maintaining a constant blood volume. Isovolumetric contraction occurs following diastole, during which the heart muscle contracts, increasing the pressure until the aortic valve opens to start systole. Isovolumetric relaxation occurs following systole, during which the heart muscle relaxes, decreasing the pressure until the pulmonary vein valve opens. In our model, which instantaneously switches from diastole to systole, the heart pressure P h is determined in a way that is consistent with the isovolumetric constraint. In particular, during systole, and during diastole, (3.10) Consequently, while P a and P v are continuous functions of time, P h experiences jump discontinuities at the transitions between diastole and systole. Simulations can be done by sequentially integrating the systolic and diastolic equations, and repeating. With model parameter values taken from the table, the simulations of the model equations exhibited realistic values. As shown in Fig. 3, the stroke volume (the amount of blood pumped out of the heart in one heartbeat) was found to be approximately 70 mL, which is consistent with typical values. The aortic pressure P a varied between approximately 70 mm Hg (diastolic) and 120 mm Hg (systolic) within a typical period of 1 second, as shown in Fig. 3. These values are within the physiological range and the trend is consistent with physiological expectations. The pulmonary venous pressure P pv , which is a surrogate for the atrial pressure, varies between approximately 2.5 mm Hg and 22.1 mm Hg. These values are also within the physiological range. The trends are reasonably consistent with expectations given our model assumptions.
Pulmonary vascular resistance R p is known to increase in Fontan failure and an increase in this resistance is known to lead to a decrease in cardiac output. This model can be used to demonstrate the impact of pulmonary vascular resistance on cardiac output. Figure 2 shows the change in average cardiac output as a function on pulmonary vascular resistance for two different heart rates. As expected, the average cardiac output decreases with increasing resistance. What is interesting is that there is a change in slope of the curves at an inflection Distance from the heart (cm) Pressure (mmHg) Figure 5: Blood pressure values (systolic (red), diastolic (green), mean (blue)) for a healthy Fontan patient as a function of distance from the heart. point corresponding to R p ≈ 3.60 mm Hg/min/L for a heart rate of 60 beats/min and R p ≈ 3.65 mm Hg/min/L for a heart rate of 120 beats/min. After this inflection point, the cardiac output decreases more quickly for a given change in resistance than before this inflection point, suggesting that something changes at this point with regards to Fontan failure and this change is dependent on heart rate. Pulmonary resistance also has an impact on the cardiac pressure-volume curve. As shown in the left panel of Fig. 4, the cardiac pressure-volume curve shifts to the left (i.e. decreased cardiac volumes) for the case of high pulmonary resistance. What this means is that the basal volume of blood in the heart has decreased as a result of this increase in resistance. For the present case, the basal volume of the heart has decreased to almost zero, implying that a further increase in resistance would result in an insufficient amount of blood returning to the heart, which consequently would reduce the cardiac output. Conversely, we see the opposite effect in the right panel Fig. 4, which shows the cardiac pressure-volume curve shifted to the right (i.e. increased cardiac volumes) for the case of high heart rate. By increasing the heart rate for a fixed stroke volume, the cardiac output would increase, resulting in an increase in the amount of blood being pumped to the body and returning to the heart. Figure 5 shows the pressure drop as a function of distance from the heart for healthy and failing Fontan patients based on clinically measured pressure catheter data. This figure illustrates that the majority of the pressure drop occurs near the heart in the systemic arteries and that furthest away from the heart, in the total cavopulmonary connection and the pulmonary circulation, the pressures are low and nearly constant.

Spatially inhomogeneous PDE models of blood pressure distribution
In this section, our ODE approach to modelling blood flow in the Fontan circulation is extended to a PDE model. The PDE model has the advantage of allowing for spatial variation of model parameters such as compliance and resistance. This will allow for more personalization of the model to an individual patient, which should improve the accuracy and predictive capabilities of the model. Furthermore, with piecewise constant values of the model parameters, the PDE model should show a similar behavior to the ODE model, giving us some assurances as to the fidelity of the PDE approach.
To model the circulatory system as a continuous flow network in a resistive compliance vessel, we assume that blood flow is a Stokes flow, i.e. the Reynolds number is sufficiently small to allow us to neglect inertial effects. Consequently, the flux in a cylindrical tube is a Poiseuille flow for which where P is the local pressure, A is the cross-sectional area, and µ is the fluid viscosity. Now we assume that a vessel is a linear compliance vessel with A = A 0 + C P , where C is the compliance. This leads to a flux relationship for a single vessel If we have a total number of N parallel vessels all with cross-section area A, the flux is Notice that in the limit C → 0, this reduces to Ohm's Law (as it must) is the resistance per unit length. When combined with the conservation law (the total volume of circulating blood is conserved) ∂A ∂t + ∂Q ∂x = 0, (4.5) this yields which is a nonlinear parabolic partial differential equation for P (x, t) with spatially variable coefficients. In general C = C(x) ≥ 0, A 0 = A 0 (x) ≥ inf(A 0 ) > 0, and N = N(x) ≥ 1.

Boundary conditions for normal circulation
For the normal circulation (see Fig. 6), there are two hearts, the left and the right hearts, each of which satisfy a volume-compliance relationship of the form The basal volumes and compliances are different during systole and diastole periods, that is (4.9) Note that in equations (4.8) and (4.9), the subscripts l and r refer to the left and right hearts, the subscripts d and s refer to diastole and systole. In this two heart model, the 7 pressures are systemic venous pressure P v , vena cava pressure P V , pulmonary arterial pressure P pa , pulmonary pressure P p , pulmonary venous pressure P pv , aortic pressure P A , and systemic arterial pressure P a . During systole the input valves (mitral and tricuspid) are closed and output valves (pulmonary and aortic) are open, while during diastole the opposite is the case. For the normal circulation, we let 0 < x < L r be the systemic circulation, and L r < x < L l be the pulmonary circulation, and of course the domain 0 < x < L l is periodic.
For convenience, we use the following notation, during systole: systemic pressure is P 1 s := P A , pulmonary pressure is P 1 p := P pa and during diastole: systemic pressure is P 2 s := P pv , pulmonary pressure is P 2 p := P V . Systemic and pulmonary pressures are everywhere continuous functions except two points x = L l = 0 and x = L r where discontinuity jumps correspond to the left and right heart pressure jumps during the switch between systole and diastole phases.

Partial differential equation:
C(x) ∂P ∂t = q(x, P ) ∂P ∂x . (4.10) Assume that at the initial time t = 0 the pressure is P (x, 0) = P 0 (x). Systemic circulation model P 1 s is given by the partial differential equation (4.10) and conditions: Problem P 1 s : 0 < t < t 1 , 0 < x < L r , Initial data : P 1 s 0 (x) = P 0 (x), Boundary conditions for P 1 s : B 1 s (t) := P 1 s (0, t), Dynamic Flux BC at x = 0 : C ls dB 1 s dt = q(0, B 1 s ) ∂P 1 s ∂x (0, t), Neumann BC at x = L r : P 1 s x (L r , t) = 0. (4.11) Pulmonary circulation model P 1 p is given by the partial differential equation (4.10) and conditions: (4.12) These boundary conditions are defined in such a way as to ensure conservation of the total volume during systole. Indeed, from (4.5), (4.8), (4.9) and the dynamical boundary conditions it follows directly that d dt V T = 0.
4.3 Initial-boundary value problem (switch from systole to diastole at: t = t 1 ).
Initial data for the systemic circulation model P 2 s (and for B 2 s (t) := P 2 s (L r , t)) at the beginning of the diastole: (4.13) Initial data for the pulmonary circulation model P 2 p (and for B 2 p (t) := P 2 p (L l , t)) at the beginning of the diastole: (4.14) 4.4 Initial-boundary value problem (diastolic regime: t 1 < t < t 2 ).
Systemic circulation model P 2 s is given by the partial differential equation (4.10) and conditions: Problem P 2 s : t 1 < t < t 2 , 0 < x < L r , Initial data : P 2 s (x, t 1 ), Boundary conditions for P 2 s : B 2 s (t) := P 2 s (L r , t), Neumann BC at x = 0 : P 2 s x 0, t) = 0, Dynamic Flux BC at x = L r : C ld dB 2 s dt = q(L r , B 2 s ) ∂P 2 s ∂x (L r , t).
(4.15) Pulmonary circulation model P 2 p is given by the partial differential equation (4.10) and conditions: Problem P 2 p : t 1 < t < t 2 , L r < x < L l , Initial data : P 1 p (x, t 1 ), Boundary conditions for P 2 p : B 2 p (t) := P 2 p (L l , t), Neumann BC at x = L r : (4.16) Again, it follows directly from (4.5), (4.8), (4.9) and the dynamical boundary conditions that during diastole d dt V T = 0.
(4.17) Initial data for the pulmonary circulation model P 1 p (and for B 1 p (t) := P 1 p (L r , t)) at the beginning of the systole: , for x ∈ (L r + ǫ, L l ) initial data, P 1 p (L r , t 2 ) = B 1 p (t 2 ), and P 1 p (x, t 2 ) is smooth for x ∈ (L r , L r + ǫ) interpolation.
(4.18) The problem is periodic in time (i.e systolic and diastolic regimes are repeated).

Boundary conditions for Fontan circulation
The Fontan blood flow circulation has only one heart, so the model has only a single loop 0 < x < L l . The basal volumes and compliances of the heart are also different during systole and diastole periods, that is For convenience we use the following notation, during systole: blood pressure is P 1 := P A and during diastole blood pressure is P 2 := P pv . 4.7 Initial-boundary value problem (systolic regime: 0 < t < t 1 ).
(4.23) The problem is periodic in time (i.e systolic and diastolic regimes are repeated). It is a direct computation to verify that with these boundary conditions d dt V T = 0.
5 Well-posedness analysis of spatially inhomogeneous PDE model 5.1 Existence and uniqueness of the nonnegative solution in normal case.
In this section we find restrictions on the parameters of the PDE model with dynamical boundary conditions for the case of the normal blood circulation for which a nonnegative solution exists and is unique. First, we show that the pressure P (x, t) stays bounded for any time t, then we obtain uniform in time bounds for the time derivative and for the gradient of P (x, t) for some range of parameter values. Second, we employ Grönwall lemma to prove uniqueness. We start by introducing the following notations , Ω r := (0, L r ), Ω ǫ r := (L r − ǫ, L r ).
Our main result establishes parameter ranges for which non-negative solutions exist, as follows.

Proof of Theorem 1
Note that the equation (5.1) becomes degenerate if P = − inf A 0 (x) C(x) . Hence, we start by constructing a sequence of positive approximations of initial data P 0n > 0. We can choose for example P 0n (x) = P 0 (x) + 1 n (if n → ∞ then P 0n (x) → P 0 (x)). These approximations allow us to apply the theoretical background developed for uniformly parabolic equations to our problem. By taking n → ∞, as a limit, we obtain a weak solution P (x, t). We omit some technical details and derive only a priori estimates which imply the existence of this weak solution. We need to specify conditions on the smoothing functions f (x, t) andf (x, t) such that total volume is conserved on the whole time interval [0, t 2 ].
The proof of Theorem 2 is similar to the one of Theorem 1. Finally, to get well-posedness for the whole interval (0, L l ), the restrictions on the parameter values obtained in Theorem 1 should be combined with the restrictions obtained in Theorem 2.

Existence and uniqueness of the nonnegative solution in Fontan
case.
The proof of Theorem 3 is similar to that of Theorem 1.

A On super-and sub-solutions
The second order parabolic equation admits super-and sub-solutions. Here we show how to constract them for a class of bounded initial values P 0 (x). By the comparison principle this implies that P sub ≤ P ≤ P sup . First of all, we consider the problem for the sub-solution: where d 2 1 = k 1 inf C(x) , k 1 = inf[q(x, P )], with initial and boundary conditions P 1 (x, 0) = inf P 0 (x) 0, (A.2) P 1,t = k 1 C ls P 1,x at x = 0, P 1,x = 0 at x = L r (A.3) for all t ∈ (0, t 1 ); P 2,x = 0 at x = 0, P 2,t = k 1 C ld P 2,x at x = L r (A.4) for all t ∈ (t 1 , t 2 ). The problem (A.1)-(A.3) has a particular solution where λ 1 satisfies tan( λ 1 d 1 L r ) = − d 1 C ls k 1 λ 1 . The problem (A.1), (A.4) has a particular solution P 2 (x, t) = c 0 cos( λ 2 d 1 x)e −λ 2 2 t , where λ 2 is the solution of the following equation tan( λ 2 d 1 L r ) = d 1 C ld k 1 λ 2 .