Abstract
This paper is concerned with the mathematical modeling of a severe and common congenital defect called hypoplastic left heart syndrome (HLHS). Surgical approaches are utilized for palliating this heart condition; however, a brain white matter injury called periventricular leukomalacia (PVL) occurs with high prevalence at or around the time of surgery, the exact cause of which is not known presently. Our main goal in this paper is to study the hemodynamic conditions under which HLHS physiology may lead to the occurrence of PVL. A lumped parameter model of the HLHS circulation has been developed integrating diffusion modeling of oxygen and carbon dioxide concentrations in order to study hemodynamic variables such as pressure, flow, and blood gas concentration. Results presented include calculations of blood pressures and flow rates in different parts of the circulation. Simulations also show changes in the ratio of pulmonary to systemic blood flow rates when the sizes of the patent ductus arteriosus and atrial septal defect are varied. These changes lead to unbalanced blood circulations and, when combined with low oxygen and carbon dioxide concentrations in arteries, result in poor oxygen delivery to the brain. We stipulate that PVL occurs as a consequence.
1. Introduction
Hypoplastic left heart syndrome (HLHS) is a congenital heart defect (CHD) in which the left side of the heart is severely underdeveloped. Without early surgical palliation HLHS is universally fatal [1]. HLHS results from a failure of the aortic or mitral valve to form. Lack of antegrade flow through the valves causes insufficient growth of both the left ventricle and the ascending aortic arch. A typical HLHS heart is compared with a normal heart in Figure 1.
(a)
(b)
Underdevelopment of the left ventricle-aorta complex resulting in critical aortic valve stenosis or aortic valve atresia with an intact ventricular septum is the most recognized form of HLHS. There are corresponding changes in the right side of the heart in the case of HLHS. All right sided cardiac structures are larger than normal including the right atrium, pulmonary artery, and pulmonary valve [2].
As blood returns from the lungs to the left atrium, it must pass through an atrial septal defect to the right side of the heart. In cases of HLHS, the right side of the heart must pump blood to the body through a patent ductus arteriosus (PDA). This maintains fetal parallel circulation where the right ventricle is the only active pump [3]. But, since the ductus arteriosus (DA) usually closes within eleven days after birth, for an HLHS baby, blood flow is severely compromised leading to very low circulation and possible death. Hence, the management of neonates with HLHS is very complex. Treatment generally commences with vigorous infusion of prostaglandin to prevent the PDA from closing. However, reduction in pulmonary resistance after birth results in an unbalanced circulation where most of the blood goes into the pulmonary circulation thereby compromising the systemic oxygen supply. A typical way to treat this disease is shunt surgery. A shunt is a surgically created connection between the systemic arterial circulation and the pulmonary arteries. Authors in [4] present a good survey of shunt modeling and its application in planning the surgery.
The need for a detailed study of the problem of HLHS stems from our recent studies [5–7] to predict the occurrence and extent of periventricular leukomalacia (PVL) after Norwood surgery in neonates. The PVL is a form of white matter brain injury, characterized by necrosis (more often coagulation) of white matter near the lateral cerebral ventricles. It can affect newborns and (less commonly) fetuses; premature infants are at the greatest risk of the disorder. Affected individuals generally exhibit motor control problems or other developmental delays, and they often develop cerebral palsy or epilepsy later in life [8]. Despite being full-term infants, PVL is found in more than 50% of the neonates after cardiac surgery [1, 9]. As there has been a dramatic reduction in mortality rates following surgery for complex CHDs, there has been an increasing recognition of adverse neurodevelopmental sequelae in some survivors. Evaluation of children following neonatal repair of CHD demonstrates long term neurodevelopmental morbidity and learning disability in survivors of these infant heart surgeries. Management strategies before, during, and after surgery including the type of support during operation have been implicated as factors in postoperative neurodevelopmental dysfunction. Deficiency in oxygen circulation and low carbon dioxide concentrations () have also been considered as important factors affecting the morbidity rates of neurodevelopmental dysfunction [6].
In our previous work in the field of PVL prediction, we have applied computational intelligence techniques to the data obtained from 59 patients which were collected at the Children’s Hospital of Philadelphia (CHOP). Computational intelligence (CI) is a set of nature-inspired computational methodologies, examples of which include Fuzzy Logic systems, Neural Networks, and Evolutionary Computation. The decision trees that we have developed suggest some ranges for critical hemodynamic parameters, such as oxygen concentration and blood pressure, which predict a high probability for the occurrence of PVL [5, 10, 11]. But since these numbers have been obtained from a retrospective review of limited set of subjects, they cannot be used as a general applicable criterion for all neonates. Hence, it is important to carry out physiological modeling of the cardiovascular system to begin to understand the underlying causal relationships between a particular set of parameters with the occurrence of PVL.
In this paper a lumped parameter model of the HLHS circulation has been developed to study blood pressure and flow in different parts of the cardiopulmonary circulation system. The diffusion modeling of oxygen (O2) and carbon dioxide (CO2) is also included in the model. Changes induced in O2 and CO2 concentrations and variation of blood flow rates in different parts of the body due to changes in some important model parameters have also been studied.
2. Pressure Flow Model
Several lumped parameter models of the cardiovascular system have been developed in past research, for instance, [12–14]. This paper differs from previously developed cardiovascular models in that this model is tailored to the unique physiology of the HLHS heart in comparison to the abundant literature which deals with the simulation of the normal heart. Furthermore, this model combines blood gas modeling with the circulation model to represent a more comprehensive model of the HLHS physiology. Moreover, by using parametric analysis in this study we have been able to track the changes in the heart for before and after surgery conditions. Like all other lumped parameter models, this model is limited by simplification assumption that it makes, for example, overlooking solid fluid interactions of arteries and blood flow and Newtonian flow assumption. For this work the lumped-parameter approach used to model the HLHS circulation is shown in Figure 2. Each box represents a lumped parameter model for a complex system of blood vessels and heart components. This model is built based on methodologies previously published for the fetal [15] and neonatal cerebral circulation [16–19] and the Norwood procedure [20, 21]. The complete model is composed of three main parts: a hypoplastic left heart, systemic circulation, and pulmonary circulation. Each compartment is made of resistances, capacitances, and inductances. Resistances are used to model the flow in the arteries and veins, and capacitors are used to represent the elasticity of these vessels. Inductors are not typically used for modeling the flow in the veins because veins do not function primarily in a contractile manner and so their inductance values are considered to be negligible in comparison with that of the arteries.
2.1. Hypoplastic Heart
The heart has been assumed to be composed of three parts: right atrium (RA), left atrium (LA), and right ventricle (RV). Here, as mentioned earlier, we have considered the most severe case of the hypoplastic heart in which the left ventricle is completely blocked (the most common form of HLHS is mitral atresia/aortic atresia; atresia implies no flow through the valve). Hence the left ventricle is not taken into consideration in the model.
The activity of the heart is modeled following the treatment in [22]. For both the atria and the ventricle, the total pressure is expressed as a nonlinear function of the volume and cycle-time:where is the pressure, is volume, is the unstressed volume, and is the maximum elasticity of the heart wall during the heart cycle. The activation function is the driving force of the heart and models the release of Ca2+ which initiates the contraction of the heart muscle [23].
Since the heart muscle’s contraction is different for systole and diastole cycles, the activity function is defined by different differential equations in the systole and diastole. During diastole, the activation function is defined by the following:where is the relaxation rate. For the systole period it is defined by the following:where is the excitation rate and is the limiting value of the activation function. The solutions to the above linear differential equations yield the following:where is the systole time and is the onset of diastole which we set equal to zero assuming the heart cycle starts from onset of diastole. This assumption is just for convenience and does not impact generalization of the model. The parameter values of , , , and are different for ventricles and atria and are shown in Table 1 [23]. The activation functions are functions of time for the ventricle and the atrium and are shown in Figure 3.
The is a function of volume of the ventricle to account for the decreasing elastance with increasing volume [20]. Considerwhere is constant; is a negative-valued constant to account for decreasing elastance of the ventricle. is zero for the atria resulting in a constant , consistent with the literature. Again, is the unstressed volume of the chamber which is the volume at zero pressure.
This is the -intercept of the tangent at the end-systolic point for the pressure-volume (-) curve. For the diastole, the ventricle and atria fill through an exponential - function reflected in the equations below.
Hence, for maximum isometric pressure, we have [22]where is normalized . The atria and ventricles are modeled with variable capacitors to account for the time-dependent relationship of pressure with volume and also viscous losses. A flow resistance has been introduced between the right and left atria to account for the defect in the walls of the heart permitting a mixing of blood flow between the atria. In the present work, this mixing has been modeled as an orifice unlike [24] where it was modeled as a simple resistor. The reason for modeling ASD as a variable nonlinear resistor instead of as a simple linear resistor is because of its very small diameter in comparison with other compartments which leads to a nonnegligible nonlinearity. A nonlinear pressure flow relationship (Darcy-Weisbach equation) is used as follows:
The tricuspid and the pulmonary valves are also modeled as unidirectional orifices (diodes) and a similar pressure flow relationship has been used for them. This model represents an idealized situation since, in reality, the tricuspid valve does leak sometimes and the flow is not completely unidirectional.
2.2. Pulmonary Circulation
The pulmonary circulation is divided into three compartments: proximal pulmonary arteries (PA), pulmonary arterial bed (PAB), and pulmonary venous bed (PVB). The PA and PAB are modeled using and PVB by circuits. All these circuit elements have constant values. Flow enters into the PA from the pulmonary valve and the blood flows out to the left atrium (LA). It should be noted that the PDA is present in the newborn which normally closes after 5–10 days. To model it, a nonlinear resistor is added between the pulmonary artery and the aorta. The addition of the PDA to the model is another major difference between this model and other previously developed models; this is especially important in the current study because of our focus on HLHS.
2.3. Systemic Circulation
The systemic circulation is divided into four compartments: descending aorta (DA), systemic arterial bed (SAB), systemic venous bed (SVB), and systemic large veins (SLV). The DA and SAB are modeled by and SVB and SLV are modeled by circuits. Blood comes in from the PDA and flows to the right atrium.
2.4. Fluid Flow
For each compartment, time-dependent variables, such as pressure in the capacitance, are expressed as differential equations and since these equations cannot be solved analytically, numerical integration is used. This approach calculates solutions with variable time steps to avoid mathematical instability that would lead to divergence of the solution. For every RLC compartment shown in Figure 4, analysis has been carried out as shown.
The relationship between pressure and volume in the capacitor leg in Figure 4 can be calculated using the following:where is the change in the compartment volume calculated from the pressure differential. The pressure change in Figure 4 is given by the following:Hence, for every compartment, the change in the output flow rate is expressed using the following equation:ASD is assumed to be a variable nonlinear resistor; therefore, we can write the relationship between and pressure change through Darcy-Weisbach equation: where and are the right and left atrial pressures, respectively. By assuming realistic initial values, numerical integration is carried out for a sufficient number of heart cycles to achieve steady-periodic values for every parameter. The PDA flow rate is estimated from the following:where and are the pulmonary artery and descending aorta pressures, respectively.
2.5. Oxygen and Carbon Dioxide Diffusion Modeling in the Vascular System
The oxygen diffusion takes place from lungs to the alveolar capillaries and then from the blood capillaries to the body tissues. At steady state, the oxygen release and uptake will be equal, and hence we need to model only one diffusion process to find the partial pressures. The process of oxygen diffusion modeling is the same as what has been done in [20] but carbon dioxide calculation is added in our model. The reason for including CO2 diffusion to the model is that our previous findings have suggested blood CO2 concentration to be an important parameter in prediction of the PVL occurrence.
On the whole, the oxygen uptake and release equations can be written as was done in [25]. Considerwhere is the oxygen uptake rate in the lungs expressed in mL/min, and it is specified as an input to the system, is the whole body oxygen consumption, is oxygen concentration in the pulmonary vein, is oxygen concentration in the aorta, and is oxygen concentration in the veins. Since the PDA concentration of oxygen in aorta and pulmonary artery is equal, we simply replace oxygen content in the pulmonary artery with its value in the aorta in the second equation. It should be noted that this assumption is only true for the MA (atresia) patients.
The term is the average pulmonary flow obtained by taking an average of (atrial septal defect flows) over multiple cardiac cycles. There are two reasons. Firstly, under steady state conditions, all the pulmonary venous return flow passes through the ASD and pulmonary venous flow is equal to pulmonary flow. Hence, on average atrial septal defect flow and pulmonary flow are equal. Secondly, a large fraction of (ideally half) the pulmonary artery flow goes to aorta, and hence, to find the average pulmonary flow, it is more reliable to calculate . Similarly, is the average systemic flow obtained by taking an average of PDA over the heart cycles.
Oxygen exchange in alveoli can be modeled using Hill’s approach [26]. Oxygen in the blood is determined by the amount of oxygenated hemoglobin (98%) and the oxygen dissolved in blood (2%). So, the total volume of blood oxygen concentration () depends on the saturated blood oxygen () and the partial pressure of blood oxygen () by the model from [16]; this is shown in (15). The term is the concentration of O2 in the blood hemoglobin at 100% saturation and the term is concentration of dissolved O2. and have units of mL·mL−1mmHg−1 and the unit of concentration is mL/mL. One hasDifferentiating (15), we get the following:
The relationship between oxygenized hemoglobin and dissolved oxygen is given by the oxygen dissociation curve: where HbO2 is the oxygenized hemoglobin, is the partial pressure of oxygen at 50% saturation (), and is a constant. Differentiating the above equation, we get the following:where By combining the above equations, we have the following:Assuming that oxygen in the capillaries of the lungs is at a constant partial pressure of 150 mmHg and that there is a continuous oxygen transfer from the lungs to the pulmonary blood capillaries, if we take the volume in the capillaries as a single unit with a volume , then the time available for the transfer of oxygen to blood capillaries can be assumed to be the time taken for the pulmonary flow to fill the capillary volume, . One haswhere is the pulmonary flow rate. At any point of time, the diffusion flux from the lungs to the blood capillaries is given by the following:where is the total lung diffusion capacity for oxygen. The diffusion flux of oxygen will lead to a change in the concentration of the volume in the blood capillaries, and hence we could writeThus, we have the following:Substituting for we getHence, at any instant, the change in the partial pressure of oxygen in blood capillaries by diffusion is given by the following:Integrating this equation over the limits, by using the numerical trapezoidal integration method with a time step of 0.01 s,where is the partial pressure of oxygen in the pulmonary veins. Since this is a complex function of , integration is carried out numerically. An estimated value of is used as an initial guess and the value of is calculated from the integral. If (27) is not satisfied, an improved value of is chosen and integration is repeated. This procedure is continued until (27) is satisfied to the required level of accuracy. The values of and are checked to satisfy the integral equation and if the left hand side of (27) turns out to be less (or more) than the right hand side, then the integration is carried out again with a lower (or higher) value. This process is repeated until the equation is satisfied to the required level of accuracy. Using the final value of and oxygen release equation, the value of is calculated.
The exchange modeling for CO2 has been implemented using the same approach as we used for oxygen. The only difference is that the CO2 is released in the lungs and it is taken up by the blood from the body tissues. Hence, we havewhere is the CO2 release rate in the lungs expressed in mL/min, and it is specified as an input to the system, is the whole body CO2 release rate, is CO2 concentration in the pulmonary vein, is CO2 concentration in the aorta, and is CO2 concentration in the veins.
Referring to [17, 18] for the relationship of concentration of CO2 to its partial pressure at a temperature of 37°C and pH = 7.4 and neglecting the effect of oxygen saturation (which is reasonable), we getwhere is the blood CO2 concentration in mL/mL. The diffusion equation for CO2 becomeswhere is diffusion capacity of lung for CO2 which is known and is alveolar pressure of CO2, which is set to 41 mmHg. Using (30) and (31), and after integration, we haveEquation (32) is numerically integrated with an initial value of such that the value of satisfies carbon dioxide release equation in the lungs. It is iterated until the equation is satisfied with a certain level of accuracy. Again, using the value of and using the equation of carbon dioxide taken up by the blood, we get a solution for .
2.6. Parameter Values
The parameter values for the case of a newborn HLHS patient were taken from a child after stage 1 circulation with a neoaorta and a Blalock-Taussig shunt replacing the PDA [20]. Since the first step Norwood procedure is applied for a child in his/her early few weeks, the parameter values for these infants would be comparable; these values are listed in Tables 1–3. It should be noted that , , and are the only variables which change postoperatively; we perform an analysis by varying the parameters to include a wide range of patients.
The value of the resistance for the pulmonary arteries is almost a hundred times that which is used in [20]. This is an estimated value to show very severe conditions. This is because, for a newborn baby, it is well known that pulmonary resistance is very high, falls precipitously after the first breath, and continues to fall over the first few weeks of his life. The average body surface area (BSA) for an adult male is around 1.5 m2 and the average surface area for a newborn child is 0.30 m2. Hence a scaling factor of five can be assumed while estimating some of the parameters for a child from the parameters of an adult. Although this method is widely used by clinicians in estimating the parameters, it brings some uncertainty to the model; hence, a parametric analysis has been carried out in order to understand the effect of those variables on the outputs.
The oxygen consumption rate in an average adult is 300 mL/min [24] and hence the average oxygen consumption rate in a child can be assumed to be approximately 60 mL/min. The total lung diffusion capacities for O2 and CO2 are 62 mL/min (mmHg)−1 and 478 mL/min (mmHg)−1, respectively, for an adult. Considering the scaling factor of 5 and the fact that the diffusion capacity is inversely proportional to the square of the surface area, we have a diffusion capacity for a child for O2 and CO2 as 2.5 mL/min (mmHg)−1 and 19 mL/min (mmHg)−1, respectively. Average CO2 release rate by an average adult is 400 mL/min, and hence the average release rate for a child is 80 mL/min.
2.7. Numerical Methods and Algorithms
The simulations and modeling presented in this paper have been done using Matlab. We used the trapezoidal integration method with a time step of 0.01 s for numerical integration. The numerical differentiation for solving ordinary differential equations presenting the HLHS model has been done using explicit fourth-order Runge-Kutta method [27]. In our simulations we have set the time step and time span to 0.01 s and 10 s, respectively. We have then ignored the first few cycles to drop the transients and to achieve steady-state solution. The initial condition has been set to 0.1 of the respective units for the pressures and flows. The computer used to run the simulations has a 32 GB of RAM and an Intel Xeon X5680 processor. The computation time to run the simulations for 10 s was 6.19 s.
3. Results
Numerical results were obtained using the parameter values as discussed above. Typical important pressure curves for various points of the body are shown in Figures 5 and 6. The pressures depicted by the model in Figures 5 and 6 are consistent with the atrial, ventricular, and pulmonary pressures of the HLHS heart. The pressure values of the different sections of the heart through the two contraction cycles are consistent with a heart that has a blockage/defect of the left ventricle as in patients with hypoplastic left heart syndrome, where the right ventricle is responsible for circulating blood to the lungs as well as throughout the body.
The variations of flow rate in the aorta and the pulmonary artery are shown in Figure 7. Another very important curve is the - loop for the right ventricle, which is shown in Figure 8. There are no discontinuities or nonsmoothness in the figure since all the relationships that are used for calculating the pressures and the volumes are exponential. Note that the heart rate in this paper is assumed to be 160 bpm which is calculated from the mean values of patient data. The end-diastolic and systolic volumes for the right ventricle are 49.9 and 18.2 cm3, respectively, with an ejection fraction of 63.5%.
Figure 9 displays the PDA flow rate. Increasing the resistance decreases the flow rate as might be expected.
3.1. Parametric Analysis
We performed a parametric analysis to investigate the effect of some clinically important parameters on the model performance and also to provide a way to compare different stages in the HLHS baby’s life. When a baby is born, the PDA is fully open; hence is low. When the baby takes his first breath, pulmonary resistance falls sharply in the first few hours but still remains at a high value for the first few days of life as it slowly drops down to normal. To model this phenomenon, we assume and mmHg·s·mL−1; the resulting pressure curves for the different parts are as shown in Figure 10.
As can be seen in Figure 10, aortic pressures are in the normal range of 75–115 mmHg and the pulmonary pressures are also normal. But as the baby ages by 5–10 days, the PDA begins to close and hence increases.
As the pulmonary resistance further drops to mmHg·s·mL−1, less flow crosses PDA to aorta. In this situation, since the pulmonary pressure is higher than the aortic pressure, there would be more systemic flow than pulmonary flow. The pressure curves for this stage are shown in Figure 11. In this situation, the condition becomes extremely dangerous. All the flow goes to the lungs and there is a reverse pressure gradient across the PDA such that the pressure in the PA is less than the pressure in the aortic arch and there is no flow going to the vital organs. Pressure builds in the RA because of the extra small pulmonary flow and death ensues.
Norwood procedure is performed to create a neoaorta from the common pulmonary artery as the cardiac outflow track. A synthetic shunt of constant diameter is placed to stabilize the pulmonary blood flow and the PDA is legated.
The constant diameter of the shunt limits the amount of flow to the low resistance PA, thereby stabilizing , a model which has been presented in [20].
Three parameters,, , and , were chosen and were varied over a range of physiological resistance, and the resulting changes in the hemodynamic parameters were analyzed. Varying ASD resistance () corresponds to varying magnitudes of the septal defect. Varying the pulmonary artery resistance () corresponds to the increasing age of the baby. Finally, the last parameter varied is the resistance within the PDA (). This parameter can be controlled by an infusion of prostaglandin to prevent PDA closure.
Two parameters are varied in each plot to study the effects of the parametric variation on systemic hemodynamics. The graphs of important hemodynamic parameters such as cardiac output, pulmonary to systemic flow ratio, and oxygen and carbon dioxide partial pressures are presented (Figures 12–15). Systemic oxygen delivery as a function of increasing is shown in Figure 12. We can observe from this plot that oxygen delivery will decrease with increase of both and . Although this is a fairly obvious fact, the interesting finding is that in the case of very high the oxygen delivery is not monotonically decreasing, and it has a maximum for some intermediate values of . We can infer that in the case of neonates with very high an increase in to some extent could be helpful.
(a)
(b)
(a)
(b)
(a)
(b)
As can be observed from Figures 13(a) and 13(b), the pulmonary-to-systemic flow ratio () changes with the change in the model parameters. As increases, this ratio decreases because the pulmonary flow keeps on decreasing; this is because, at steady state, all the pulmonary flow passes through the ASD, and an increase of leads to a decrease in the flow. As increases, the systemic flow decreases since all the systemic flow passes through ductus arteriosus which leads to an increase in the ratio.
A ratio around one is optimal for the body [28] and hence a suitable area for the flow can be selected from these graphs accordingly. Considering the parameters that have to be changed, such a point can be reached. The direction of blood flow depends upon the resistance, and hence, in the case of HLHS, blood follows the path of least resistance. In high resistances in constant , the ratio decreases because, due to high pulmonary resistance, the amount of pulmonary flow decreases and systemic flow increases leading to a decrease in the ratio. Plots show that in some combination of PA, ASD, and PDA resistances the optimal can never be achieved. For example, if mmHg·s·mL−1, the optimal ratio is unreachable.
As shown in Figures 14(a) and 14(b), as increases, the arterial O2 saturation also increases and becomes nearly constant at high values of the resistance; however, it decreases monotonously with . The decrease with becomes more profound at lower values of . This trend is also observed with . For constant , systemic arterial O2 concentration increases with increasing .
As increases, cardiac output (sum of pulmonary and systemic flows) decreases slowly as the systemic flow keeps on decreasing and this counters the effect of the increasing pulmonary flow as shown in Figures 15(a) and 15(b). As increases, cardiac output decreases significantly due to a rapid reduction in the pulmonary flow. Hence, it can be deduced that the sensitivity of the pulmonary flow to is very high. A similar trend can be observed with playing the same role as that of .
Figure 16 shows a plot of (partial pressure of carbon dioxide) for constant mmHg·s·mL−1 and varying values of and . We focus on since our previous studies based on computational intelligence (CI) techniques have shown that 43 mmHg is an important threshold for predicting PVL. By inclusion of in our model we attempt to find a reason for this interesting but somewhat inexplicable result of the CI technique. Our results show that all values of are less than 43 mmHg for mmHg·s·mL−1 and mmHg·s·mL−1. Considering the fact that to some extent is measureable and controllable post- and preoperatively, this finding could well lead to a valuable insight.
4. Conclusion
The goal of this paper has been to develop a lumped parameter model of HLHS, a fairly common congenital heart disease. The need for this physical modeling comes from our ultimate goal to predict and prevent one of the consequences of this abnormal circulation, namely, a form of brain injury called periventricular leukomalacia (PVL).
The exact causes of PVL are not well understood and in this paper our focus has been on understanding how some critical parameters such as flow resistances at the atrial septal defect or patent ductus arteriosus might alter systemic flow and systemic oxygen delivery, which could be reasonably expected to have a causal relationship to PVL occurrence.
In this paper, we have compared the results of HLHS circulation (preoperative condition) to Norwood circulation (postoperative condition). Results show the manner and extent of changes in the pressure and blood flow in different parts of the body. In reporting the modeling results, we have mostly focused on the pulmonary to systemic flow ratio (), O2 delivery, and O2 saturation because we believe that these parameters are important risk parameters of the PVL occurrence. By use of the model, the following points were demonstrated or confirmed.
Different sizes of the PDA and the ASD and different resistances of the PDA and the ASD will result in changing the regulation of pulmonary and systemic flow and hence these resistances play a critical part in our goal to predict the occurrence of the PVL.
The direction of blood flow depends upon the PA resistance and that of the PDA resistance. Increase in the PDA resistance leads to blood flowing where the resistance is the least. This results in an increased ratio. With an increase in the PDA resistance a greater portion of the cardiac output goes into the lungs. This will lead to a systemic hypoperfusion and poor O2 delivery. This result is consistent with the results presented in [29].
By comparing the slope in cardiac output and pulmonary to systemic flow variation plots, one could observe that the rates of change (slopes) are higher for variation of PA resistance in comparison with the variation of ASD resistance. This means that the PA resistance is more influential on the CO and pulmonary to systemic flow ratio than the ASD resistance. This suggests that the PA resistance is more directly involved in the determination of the ratio and the cardiac output.
The optimal O2 delivery is achieved when balanced pulmonary and systemic perfusion is established, namely, when is around 1. Our results show that, for some values of and , the optimal can never be achieved. Since insufficient O2 delivery is considered to be a potential cause for PVL, it seems reasonable to infer that these values of and lead to high risks of PVL.
It is known that the correlation between venous O2 saturation and O2 delivery is better than the correlation between arterial O2 saturation and O2 delivery. However, since measurement of the mixed venous saturation in clinical practice remains difficult, we just considered the case of arterial oxygen saturation [30] in this study.
The developed model also has some other interesting applications. For example, with some minor modifications to cover the partly underdeveloped left heart, the proposed model could form a good foundation for the analysis of newborns with underdeveloped left ventricles and model based design of left ventricular assist devices.
Given that the exact causes of PVL are still not clearly understood, the main goal of this paper has been to develop a sound mathematical model to investigate HLHS, a severe congenital deformation which has the highest correlation with the occurrence of PVL as suggested by preliminary postoperative data collected from Children’s Hospital of Philadelphia (CHOP). Also, the paper aims to provide a tool to better understand the main factors that drive the HLHS physiology by a suitable variation of physiological parameters.
Our previous computational intelligence based findings [5, 7, 10] suggest some important factors in the occurrence of PVL; however, since the exact physiological reasons and causes of those findings are not clear, the current model has been developed to discover cases that will lead to the PVL situation and to analyze them. We hope that the developed model could open the door for further investigation of the HLHS syndrome and also to its connection with PVL. The postoperative results are consistent with what is reported in [20] which makes this model potentially valuable for clinical use.
Some of the limitations of the current work are as follows.(i)The effect of respiration and regulatory mechanisms such as HR baroreflex and total peripheral resistance baroreflex on the cardiopulmonary performance is not considered in building the model.(ii)The values of the components of the analog electrical circuit are estimated and there is still no direct, reliable, and safe way to measure them with high confidence.(iii)The results of the model should be validated clinically. A parameter adaptation algorithm for patient specific modeling should be developed and implemented on the patient specific data. The correlation between the PVL occurrence and different values for ASD, PA, and PDA resistances should be calculated for the aim of PVL occurrence prediction.
To improve the model further and to add more verisimilitude with reality, it is important to include other factors into the system such as autoregulation and to investigate its effects. Our current work is continuing on this important problem and will focus on addressing some of the above limitations.
Finally we would like to comment about the general utility of the procedures developed in this paper. With all the obvious limitations, the method of constructing a mathematical model pursued in this paper is still appealing because it permits the change of parameters very easily to understand their influence on the usual key clinical outputs such as . Such models exist for other conditions but none exist for PVL based on preoperative conditions. By use of a parametric analysis on computer models we start to acquire the very useful ability to predict pre- and postoperative conditions for specific patients and also to develop patient specific models to investigate physiological states of individual patients. That said, we would like to reiterate the limitations of the current model and hope to develop and inspire development of improved models with appropriate corroboration.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The research reported in this paper is supported by a grant from National Institute of Health (no. 1 R01 NS 72338 01A1).