With circulatory pathology, patient-specific simulation of hemodynamics is required to minimize invasiveness for diagnosis, treatment planning, and followup. We investigated the advantages of a smart combination of often already known hemodynamic principles. The CircAdapt model was designed to simulate beat-to-beat dynamics of the four-chamber heart with systemic and pulmonary circulation while incorporating a realistic relation between pressure-volume load and tissue mechanics and adaptation of tissues to mechanical load. Adaptation was modeled by rules, where a locally sensed signal results in a local action of the tissue. The applied rules were as follows: For blood vessel walls, 1) flow shear stress dilates the wall and 2) tensile stress thickens the wall; for myocardial tissue, 3) strain dilates the wall material, 4) larger maximum sarcomere length increases contractility, and 5) contractility increases wall mass. The circulation was composed of active and passive compliances and inertias. A realistic circulation developed by self-structuring through adaptation provided mean levels of systemic pressure and flow. Ability to simulate a wide variety of patient-specific circumstances was demonstrated by application of the same adaptation rules to the conditions of fetal circulation followed by a switch to the newborn circulation around birth. It was concluded that a few adaptation rules, directed to normalize mechanical load of the tissue, were sufficient to develop and maintain a realistic circulation automatically. Adaptation rules appear to be the key to reduce dramatically the number of input parameters for simulating circulation dynamics. The model may be used to simulate circulation pathology and to predict effects of treatment.
mathematical models of physiological processes have become increasingly powerful because of increases in physiological knowledge and speed of computers. Consequently, in the near future, computer models may improve diagnosis and aid in prediction of therapeutic effects in patients. Quantitative diagnosis of circulatory pathology requires a model allowing patient-specific simulation of hemodynamics in heart and circulation.
Several models of the whole closed-loop circulation (20, 24, 32, 35, 43) were designed by integration of a selection of subsystems such as heart chambers, valves, and blood vessels. The large number of needed, but often unknown, parameters appears to be a major problem. Parameter values almost certainly differ between healthy individuals, let alone between patients. Consequently, clinical applications of patient-specific models are practically nonexistent. At least part of the variability is likely related to adaptation of the body to the physiological environment. Therefore, we hypothesized that the use of adaptation rules might reduce the number of parameters and facilitate patient-specific computer simulation.
During the last decade, it has become clear that, within the circulation, mechanical signaling is important. From the very beginning of its development, hemodynamic load of the heart is a key requirement for normal cardiac development (36). Hypertrophy is a known response to control mechanical load of the cardiac tissue (28). In blood vessels, shear stress along the endothelial cells regulates lumen diameter (7, 16), and wall stress is likely to control wall thickness in some way (6).
The aim of the present study was to investigate the added value that may be obtained from knowledge of circulation dynamics by a smart combination of (partly) known physiological principles in a comprehensive model of the total circulation. Therefore, we examined the extent to which modeling of circulation hemodynamics could be made more realistic by incorporation of the already known adaptive behavior of heart and blood vessels to mechanical load imposed by the circulation. We hypothesized that such a model could predict stable self-structuring of the normal heart and circulation. Furthermore, it was expected that an isolated pathological disorder, e.g., a valve defect, might initiate a sequence of adaptive responses, resulting in a combination of events that is typical for that pathology. In diagnosis of circulatory pathology, flow and pressure disorders are mutually related. By combining different measurements, patient-specific modeling may help us find the underlying disorder of the system, which may not be recognized from analysis of the separate measurements. Given a diagnosis of pathology, different treatments may be simulated to find the best one. During treatment, the specific patient can be simulated and evaluated on the basis of proper additional, preferentially noninvasive measurements. Furthermore, treatment can be optimized during therapy.
We designed the CircAdapt model to simulate the dynamics of heart and circulation, whereas geometry of the constituents was determined by adaptation to mechanical load. Four types of modules were used: compliant blood vessels, actively contracting chambers, valves with inertia, and peripheral resistances. In blood vessels and chambers, cavity volume and wall thickness changed in response to calculated mechanical load. Responses of the circulation to interventions were assessed by simulation of exercise and variation of preload. Flexibility of the CircAdapt model was shown by simulation of the development of fetal circulation 3 mo before birth, the intervention of birth, and the development of the newborn circulation until 3 mo after birth. The latter simulation was repeated, but with a tricuspid valve that was stenotic, starting 3 mo before birth. Finally, for the normal adult circulation, a sensitivity analysis of geometric and functional parameters to chronic changes in various physiological parameters was performed.
General design of the CircAdapt model.
The CircAdapt model is configured as a network consisting of four types of modules (Fig. 1) : chamber, tube, valve, and resistance. The heart contains four chambers. Because of the typical anatomy of the normal ventricles, left-to-right ventricular interaction has been simulated by a large outer chamber encapsulating a high-pressure-generating inner chamber, representing the left ventricle (1). The space between the outer and the inner chamber represents the right ventricular cavity. The large arteries and veins connected to the heart are simulated by nonlinearly compliant tubes, having nonlinear characteristic impedance for pressure-flow waves. Chambers and tubes are connected by valve elements having inertia and Bernoulli pressure losses (10). The valve element can also handle flow characteristics of orifices, such as venous outflow into atria and septal defects of the heart.
The core of the model is a system of differential equations in state variables with respect to time, similar to earlier models of the whole circulation (17, 20, 24, 32, 35). Examples of state variables are volumes of chambers and tubes and flow through valves. The model generates traces of hemodynamic variables, e.g., pressures and flows, as a function of time within the cardiac cycle.
Adaptation rules simulate macroscopic consequences of processes where mechanical signals induce gene expressions, causing the cell to change structure and/or mechanical behavior within reach of the cell. The generally complicated, intertwined mechanisms of adaptation are condensed in relatively simple rules. Cavity and wall geometry of heart and vessels adapts to normalize load in the constituting tissues.
Design of specific modules.
Modules are subdivided in two classes. The tube and chamber module are compliances, in which pressure depends on volume. Valve and resistance modules are flow ducts. In a network representing the cardiovascular system, compliances are nodes connected by flow ducts.
A tube module (Fig. 2A) represents the entrance of a compliant blood vessel capable of propagating a pressure-flow wave component added to a constant flow. The behavior of the entrance is similar to that of a compliance in series with the wave impedance (41). The wall of the compliance consists of fibrous tissue. When the model of a cavity with a fibrous wall is used (2), average extension (λ) of the fibers in the wall depends on the ratio of cavity volume to wall volume (Vcav/Vwall) as follows (1) Note that λ = 1 refers to zero cavity volume. Cavity pressure (pcav) depends on λ (2) (2) For the constitutive equations relating mean Cauchy fiber stress (σ) to λ, we used (3) At the normal physiological reference state, i.e., the mean state of deformation at rest, Vcav,ref occurs at a physiological fiber extension λref with Cauchy fiber stress σref and cavity pressure pcav,ref. Parameter k determines nonlinearity of the stiffness of vessel wall material.
The characteristic wave impedance (Z) of a tube is a resistance depending on the ratio of inertia (Ln) to compliance (Cn) by (4) For Ln and Cn (5) Parameters ρ, Acav, and lcav represent blood density, lumen cross-sectional area, and length of the vessel segment, respectively. Equations 1–5 are used to express pcav and Z as a function of Vcav (6)
A chamber module (Fig. 2A) represents a cavity with a wall that contains myofibers. Conventionally, a chamber has been described by the time-varying elastance model (26). More detailed finite-element models coupled sarcomere mechanics (14) to cavity mechanics (4, 15). From these complicated models, a simple one-fiber model was derived in which cavity mechanics were coupled to a representative sarcomere length and myofiber stress (2), similar to the above-mentioned tube element. The chamber differs from the tube, in that the relation between λ and σ depends on time, thus allowing active generation of external work.
Sarcomere length (Ls) in the wall depends on λ (Eq. 1) (7) with Ls,V(0) representing Ls at (extrapolated) zero Vcav. Sarcomere mechanics are represented by a passive elastic element (subscript “pas”) in parallel with an elastic element in series with the contractile element (subscript “sc”) with length Lsc. σ and velocity of contractile element shortening depend on time (t), Ls, Lsc, and contractility (Cs) (8) Many models of sarcomere mechanics can be written in the format of Eq. 8, with which it is possible to extend the number of state variables in the description of sarcomere mechanics. We have used (9) The function f1(t,Lsc,Cs) represents stiffness of the series elastic element. Active stress is proportional to the length of the series elastic element (Lse = Ls − Lsc). Parameters νmax and Lse,iso represent velocity of sarcomere shortening at zero load and Lse for isometric contraction, respectively. The applied contraction model is described in more detail in the appendix. Equations 1 and 2 have been used to calculate pcav from σ.
A valve module (Fig. 2B) is a flow duct where cross-sectional area depends on flow (q) and transvalvular pressure (Δp). Total Δp is the sum of the effects of inertia due to acceleration in time and the Bernoulli effect (10). Consider a valve with cross-sectional areas Aprox, Avalve, and Adist, proximal, within, and distal to the valve, respectively. Pressure drop (Δp) over the valve depends on blood velocities νprox, νvalve, and νdist by (10) where lvalve indicates the length of the channel with inertia. Cross-sectional area (Avalve) depends on q and Δp by the following algorithm (11) Aopen and Aleak are given valve cross-sectional areas of the open and closed valve, respectively. With Eq. 11, the closed valve immediately opens by a forward pressure drop. Rewriting of Eq. 10 with the use of Eq. 11 results in an expression for the time derivative of flow as a state variable (12)
A resistance (Fig. 2C), generally representing the microcirculation, connects an arterial to a venous tube module. Δp between the arterial and venous compliance determines q by (13) where Rp represents peripheral resistance. The pressure-controlled flow source (q) connects the resistive entrances of the tubes. Thus effects of tube wave impedance (Z) and Rp on flow are mathematically separated, while properties of the three-element windkessel model (41) are maintained.
Additional characteristics of the system.
First, the left-to-right ventricular interaction is described by setting up the cardiac ventricles (Fig. 1) as an outer wall encapsulating the right ventricular cavity (rv) and left ventricle (lv). The outer wall generates right ventricular pressure with the summed volume change of left and right ventricles (1). The inner wall generates the transseptal increment from right to left ventricular pressure and a volume change of the left ventricle only
Second, the atrioventricular delay of activation was 15% of cardiac cycle time.
Third, to simulate blood pressure control, circulating blood volume was adapted until mean arterial pressure reached the desired mean blood pressure level.
To simulate structural responses to chronic load, adaptation rules were implemented as a switch. If an adaptation rule was switched on, after a cardiac cycle the signals for adaptation were calculated, resulting in structural responses according to the rules as proposed below. In simulations with acute interventions, adaptation rules were switched off, thus maintaining short-term temporal invariance of structural parameters. In the steady state of adaptation, the rules below are satisfied. Per adaptation mechanism, the time constant could be adjusted with the applied adaptation fraction per adaptation cycle. In the present study, we paid more attention to the steady-state solution than to the dynamics of adaptation.
Adaptation of the tube.
In a blood vessel, Acav,ref adapts to maintain a level of flow shear rate (Se) along the endothelium of the blood vessel. With the assumption of a parabolic profile with mean flow qref (15) The vessel wall should be thick enough to withstand a maximum pressure (pmax) with maximum wall stress (σf,max), which is a property of the vessel wall material. With use of Eqs. 1 and 2 for maximum pressure load with Acav/Awall = Vcav/Vwall (16) Equations 15 and 16 hold for different conditions, i.e., mean physiological working pressure (pref) and maximum occurring pressure (pmax), respectively. If Eq. 6 is applied to both loading conditions (17) Solving the system Eqs. 16 and 17 by eliminating Acav,max for Awall (18) pmax may appear very low for veins, with Eq. 18 leading to very thin vessel walls. One should, however, realize that a living body is subject to arbitrary movements, causing pressure fluctuations in blood vessels. We assume that the body is regularly subject to an impact with velocity νimpact, e.g., due to jumping. A column of blood moving with this velocity causes a flow shock wave, resulting in a pressure step by the wave impedance Z (Eq. 6). Thus we used (19)
Adaptation of chamber geometry.
Adaptation of the chamber has been derived from a previous study on adaptation of the cardiac structure to mechanical load (3). Cs (Eq. 8) increases with maximum sarcomere length (Ls,max) during the cardiac cycle (20) For Ls,max = Ls,max,adapt, contractility is set to Cs = 1, representing the state of adaptation during exercise. Cs can slightly exceed unity by a not critical fraction a, thus enabling control by adaptation around unity. For atrial and ventricular myocardium, parameter b is 30.0 and 4.0 μm−1 respectively, showing that the atrium is more sensitive to maximum stretch.
Vwall, quantifying hypertrophy, increases with a chronically high Cs (Eq. 20) (21) Large myofiber strain during ejection dilates the cavity due to strain softening (9), represented by a decrease of Ls,V(0) (22)
Adaptation of pulmonary peripheral resistance.
Mean pressure drop (ppulm) over the pulmonary system was kept constant by adaptation of peripheral pulmonary resistance (Rp,pulm) after each cardiac cycle.
Solving the system of differential equations.
The status of the model is defined by a database-like structure, incorporating all modules with their parameters and variables. In solving the related set of differential equations, the derivatives of all state variables were described as a function of the state variables. Therefore, first, in the chambers and tubes, pressures and wave impedances (Eqs. 1, 2, 6, and 8) were calculated from the actual volumes. From the flows in the flow ducts, the time derivatives of volumes in tubes and chambers were calculated. Next, the pressures around the valves and then the time derivatives of valve flows were calculated (Eq. 12). Finally, time derivatives of state variables inside the nodes and flow ducts were calculated (Eq. 8). The set of differential equations has been solved numerically with MatLab 6.5.0 (MathWorks, Natick, MA) software using the function ODE23 with ∼400 time steps per heart cycle. After each heart cycle, the relevant adaptation rules were applied.
Initialization of heart, blood vessel, and valve geometry.
Symbols, values, units, and description of the major input parameters are provided in Table 1. The simulation starts from a set of not very critical initial estimates. For all left and right arterial and venous tubes, mean and maximum blood pressure levels were initially estimated on the basis of physiological knowledge. Wall and cross-sectional areas of the tubes were calculated with adaptation rules (Eqs. 15–19). Tube lengths (lcav), determining the value of the windkessel compliances (Eq. 5), were scaled to tube diameter in accordance with anatomic findings (Table 1). Systemic and pulmonary peripheral resistances followed from preset pressure differences and flow (Eq. 13). The cross-sectional areas of the arterial and venous valves were set equal to the cross-sectional areas of the connected tubes. Mitral and tricuspid valve cross-sectional areas were scaled to aortic and pulmonary valve cross-sectional areas, respectively (Table 1). For the venous atrial inlets, valve leak area was set equal to open area, thus converting a valve to an open orifice.
The first cardiac beats were simulated until a steady state was reached. Vascular cross-sectional area (Eq. 15) and mean working pressures (pref, Eq. 17) were adapted for the resting condition, which is considered the most common state. Next, a state of exercise was simulated by increasing systemic flow and heart rate by factors of 4.0 and 2.43, respectively. Vessel wall thickness and chamber geometries were adapted continuously for this state using Eqs. 18, 21, and 22. After ∼100–200 beats, a steady state was reached. Next, flow and heart rate were returned to the resting condition. The adaptation protocol was repeated until relative changes in mean valve flows per beat were <10−3.
Simulations after interventions.
Clinical data on human pulmonary venous flow (31) were compared with simulations in which preload changed as a result of an increase of mean systemic blood flow from 100 to 200 ml/s, leaving heart rate unaffected. To show that the model can simulate conditions of very different patients, the circulation of a fetus has been simulated from 3 mo before to 3 mo after birth. Fetal circulation (Table 2) was characterized by 1) an open foramen ovale between both atria, 2) an open ductus arteriosus between the aorta and the pulmonary artery, 3) a large pulmonary vascular resistance, 4) symmetrical left-to-right ventricle interaction (modified Eq. 14), and 5) a low systemic peripheral resistance due to inclusion of umbilical circulation and a relatively low arterial oxygen saturation. The acute changes around the time of birth were simulated by 1) lowering peripheral resistance of the pulmonary circulation and 2) elevating systemic peripheral resistance by excluding the umbilical circulation. Related parameter values (23) are listed in Table 2. In the fetal and newborn stages, blood pressure increased linearly from 5.0 to 8.5 kPa and from 8.5 to 10 kPa, respectively. Systemic blood flow increased from 8 to 25 ml/s and from 8.8 to 14.1 ml/s in the fetal and newborn stages, respectively. In the newborn stage, the pulmonary pressure declined from 4 to 2 kPa with a decay time of 1 wk. All blood vessel properties and sizes of the heart chambers were found by application of the adaptation rules used for the healthy adult. Because there is no steady state during growth, adaptation followed the mechanical stimuli with a time constant of 1 wk. The same simulation protocol was repeated, but with the tricuspid valve cross-sectional area fixed at 50 mm2, thus simulating a persisting restriction in the tricuspid flow channel during fetal and newborn development.
Finally, for the normal adult circulation, we studied sensitivity of adaptation to a change in the following parameters: mean systemic pressure (pref), maximum systemic flow during exercise (qexc), vessel wall shear rate (Se), maximum vessel wall stress (σf,max), exponent for vessel wall stiffness (k), isometric peak active stress (σf,act,max), and peak passive stress (σf,pas,max) in atrium and ventricles separately.
In Table 3, parameters on hemodynamics are shown as calculated for an adult in the steady state of adaptation. In Fig. 3, systemic and pulmonary hemodynamics are shown for cardiac beats at rest and during exercise. At rest, almost all time courses of pressures and volumes were in qualitative and quantitative agreement with physiological data. Pressure gradients across the arterial valves were forward during the first part of ejection and backward near the end of ejection. During atrial systole, pressure gradients across mitral and tricuspid valves were forward, resulting in short flow peaks. At the end of ventricular relaxation, atrial pressures exceeded ventricular pressures, causing fast early diastolic filling of both ventricles. The ratio of peak atrioventricular flow during early filling (E) to atrial contraction (A) (E/A) exceeded unity. With increase of exercise, E and A peaks were mutually approaching and, finally, fusing. In Fig. 4, pressure-volume plots and myofiber stress-strain plots of the ventricles and atria are shown for rest and exercise. For exercise, sarcomere length at the beginning and end of ejection reached the preset values (Table 1), thus showing that adaptation was complete. With exercise, all cavities dilated and generated more pressure.
In Fig. 5, clinically obtained hemodynamics (31) are compared with simulations, while preload is increased acutely (no adaptation). Dynamic changes in left atrial pressure and pulmonary venous flow agree quite well.
In Fig. 6, simulated results on hemodynamics of a fetus just before birth are presented. Left and right ventricular pressures were similar. Flow and volume changes were higher on the right than on the left side. Flow through the ductus arteriosus and the foramen ovale were negative, i.e., from right to left. Calculated parameters on geometry of the various circulatory elements are shown in Table 3. Immediately after birth, in the newborn (Fig. 6), pulmonary arterial pressure decreased with decreasing pulmonary peripheral resistance, causing reversal of ductus flow. Flow through the foramen ovale decreased considerably as a result of the unidirectional valve function of the foramen ovale. Right ventricular preload and afterload decreased. In Fig. 7, simulated pulmonary venous flow velocity is plotted as a function of time for the fetus just before birth and in the newborn at birth and at 1 wk and 1 mo after birth. Doppler measurements (13), redrawn in the same format, are shown for comparison.
In Fig. 8, parameters on hemodynamics, wall geometry, and maximum sarcomere length are shown as a function of time from 3 mo before to 3 mo after birth. In the fetal stage, flows and wall volumes increased with growth. Flow through the right heart exceeded that through the left heart. After the acute changes with birth, ductus flow became zero, while foramen ovale flow steadily decreased to zero. Wall volume of the right ventricle decreased below that of the left ventricle. Maximum sarcomere length, reflecting tissue preload, decreased during a short period after birth but soon became normal as a result of tissue adaptation. The same simulation, but with flow restriction in the tricuspid valve, is also shown. In the fetal stage, flow through the right heart steadily became lower than that through the left heart. Foramen ovale flow largely exceeded ductus flow, which remained constantly low. Problems became clear with birth. Shunt flow through the foramen ovale remained high, while pulmonary artery flow remained significantly lower than aortic flow. Simulated surgical closure of the foramen ovale restored flow balance between the right and left heart but caused a significant rise of systemic venous pressure. This pressure rise increased with time of operation after birth; i.e., closure at 2, 6, and 12 wk resulted in venous pressures of 0.22, 0.28, and 0.35 kPa, respectively. Without surgery, in the simulation, the latter pressure did not exceed 0.01 kPa.
Table 4 shows relative sensitivities of many of the variables listed in Table 3 to changes in a selection of parameters listed in Table 1. Sensitivity of quantity y to parameter x has been quantified as (Δy/y)/(Δx/x) × 100%, with reference to the reference state. Increase of contractile stress of both ventricles resulted in a decrease of wall volume in all four chambers. Increase of the maximum passive stress level in the ventricles resulted in an increase of mean venous pressures and thickness of venous and atrial walls. Increase of this stress in the atria had a comparable effect on venous wall thickness but not on the other quantities. Increase of mean systemic blood pressure (pref) caused wall thickening of chambers and blood vessels on the left side. Increase of shear rate along the vascular endothelium resulted in decrease of vessel diameter and increase of wave impedance. Increase of maximum blood vessel wall stress predominantly caused a decrease of vessel wall thickness. Increase of k (Eq. 3) increased resting venous pressures, decreased venous wave impedance, and increased arterial wave impedance. Increase of maximum systemic flow during exercise (qexc), possibly related to training effects by physical exercise, caused hypertrophy of all chambers and increase of systemic venous pressure at rest. Strikingly, E/A, defined as the ratio of peak mitral flow during early filling to that during atrial contraction, increased most selectively with qexc and, to a lesser extent, with atrial passive stiffness. E/A is often used as a clinical measure of left ventricular filling quality.
The computing time needed for a simulation was ∼7 s per cardiac beat on a personal computer. Complete steady-state adaptation of wall mass required ∼50–100 cycles.
The use of a few quite well-known adaptation rules for heart and blood vessels resulted in stable self-structuring of the circulation, rendering surprisingly realistic simulations of flow and pressure traces for left and right heart and large blood vessels. The model is unique in having so few independent parameters (Table 1). The majority of the geometric parameters were calculated by the model itself (Table 2). The modular setup is convenient in simulating aberrations in circulation, as shown by simulations of the circulation of the fetus and birth, being the acute transition from the adapted fetal state to the not yet adapted newborn circulation. Adaptation was simulated as a function of time after birth, reaching a new steady state within a few weeks.
Most pathologies of the circulation are caused by some isolated disorder, whereas elsewhere the adaptation rules remain unaffected. The CircAdapt model accommodates this behavior by design, which makes it very promising in simulating pathologies, the related long-term prognosis, and different proposed treatments. For instance, long-term effects of tricuspid valve narrowing in the fetus (Fig. 8) could be predicted realistically, together with the effect of treatment by surgical closure of the foramen ovale after birth. Thus a valuable tool may be developed to assess and optimize personal treatment.
Although a modular model governed by a set of differential equations is quite conventional, the CircAdapt model has several unique features. In calculating cavity pressure from volume (Fig. 2A), macroscopic hemodynamic load was quantitatively linked to stress and strain of the tissue in the wall. The relation between stress and strain in the anisotropic tissue was obtained from physiological experiments on isolated muscle preparations (34). For a chamber, the assumed relation is realistic (2). For a blood vessel, the assumption of extreme anisotropy of the wall tissue is also more realistic than that of isotropy. When conventional mechanics theory is applied to a tube of isotropic material, its length is known to be invariant during pressurization. However, an axially unloaded tube with an anisotropic fibrous wall lengthens with increasing diameter, similar to the lengthening reported for blood vessels (40). The stress-strain relation along the fibers is incorporated as a module, also providing the freedom to add active properties to the wall material. Thus vasomotor tone may be implemented. Other specific features of the model are the coupling of wave impedance with compliance in the nonlinearly elastic wall of blood vessels and backpropagation of pulse waves in veins.
The CircAdapt model is unique in having combined adaptation of heart and blood vessels, resulting in self-structuring of the circulation as a system. Adaptation often results from complex, intertwined cellular mechanisms (33), where sensed signals result in gene expressions, causing local responses, which affect local structure. Although the model can handle dynamics of adaptation (Fig. 8), it is mostly focused on finding the steady-state solution after a sufficiently long period of adaptation. The adaptation rules on normalization of mechanical load of heart and blood vessel wall are quite well documented; however, that compensation is often not complete. Imaging techniques are appropriate to follow geometric changes noninvasively. Because the model can estimate mechanical load of the tissue, the model may be used to determine mechanoadaptive characteristics noninvasively in experimental studies and in patients. For instance, the effect of aging (11) appears similar to the simulated effect of adaptation to a 42% increase of vascular stiffness exponent k (Table 4). In measurements, the related 35% increase of aortic pulse wave velocity, which is proportional to aortic characteristic impedance, was associated with 8% increase of left ventricular wall thickness. This finding closely agrees with our sensitivity analysis (Table 4), predicting an increase of 10%.
The simulations of mitral hemodynamics are quite realistic. As shown in Fig. 3, at rest, mitral flow has two peaks. When heart rate is elevated with exercise, the flow peaks mutually approach, and finally merge into one peak. It has been claimed that the physiological time course of mitral flow requires left ventricular suction (42). In the simulation, pressure was prevented from being negative. Nevertheless, a physiological time course of mitral flow developed, showing no absolute need for suction to explain mitral hemodynamics.
In agreement with real mitral flow characteristics, the early filling peak is higher than the peak related to atrial contraction (E/A > 1). Clinically, E/A is used as an index for quality of diastolic filling. In many cardiac pathologies, E/A decreases. According to the model, E/A is mainly determined by stroke volume reserve, i.e., the ratio of stroke volume with exercise to stroke volume at rest (Table 4). This suggestion agrees with a decrease of E/A with many types of diminished heart function, where stroke volume at rest is already close to its maximum. In agreement with this idea, E/A was reported to be above normal with increased stroke volume reserve in young highly trained swimmers (38).
The time courses of atrial pressure and pulmonary venous flow appear realistic (Fig. 5). Pulmonary venous flow velocity was backward during atrial contraction. There were two maxima: one during ventricular systole and one early in diastole. With increase of preload, atrial peak systolic backflow increased. In addition, the minimum flow peak during ventricular isovolumic relaxation became more prominent, albeit in the experiment more clearly than in the simulation. During most of the diastolic phase, flow velocity declined about linearly. Mean flow was larger in the simulation than in the experiment, because in the experiment the investigated vein contained only part of the total pulmonary effluent. Left atrial pressure showed two peaks: one during atrial contraction and one during ventricular relaxation. With increase of preload, mean pressure and pressure fluctuations increased, although the increase of pressure peak with atrial contraction was more pronounced in the simulation than in the experiment. Thus it was concluded that the dynamics of pulmonary venous flow were quite well understood. In humans, pulmonary venous flow has been followed during development from fetus to 1 mo after birth (13), as also shown in Fig. 7. It is interesting to study the effects of adaptation around birth, because immediately after birth the circulation is out of balance because of the sudden changes in hemodynamic load. Within a few weeks, the circulation adapts to the new situation. Effects of adaptation are reflected in hemodynamics. Before birth, simulation and experiment agreed on a relatively pronounced early diastolic velocity peak, followed by a decline until atrial contraction. In the simulation, dynamic phenomena around atrial contraction appeared exaggerated. Immediately after birth, mean velocity increased. Atrial contraction was marked by a steep minimum. The early diastolic peak was diminished. After 1 wk, mean flow velocity was lowered while the dynamic changes remained. After 1 mo, pulsatile changes had large amplitude relative to average flow velocity, causing backflow during atrial systole. All these phenomena also were observed in the measurements. Thus it was concluded that dynamic changes as found in pulmonary venous flow dynamics during the first weeks after birth could be explained by adaptation of heart and blood vessels to the newborn state of hemodynamic load.
The CircAdapt model simulates changes in geometry due to adaptation, as shown for the fetal and newborn circulation in Fig. 8. Interestingly, the consequences of a relatively simple pathological disorder, such as a flow restriction in the tricuspid valve tract from 3 mo before birth, can be predicted quite in detail, finally resulting in pulmonary atresia (37), i.e., retarded development of the right heart. One may suggest that a primary stenosis in the pulmonary valve would have similar consequences, but that is not the case. Our simulation showed that, in the case of pulmonary stenosis, the right ventricle would adapt to generate a very high pressure with a thick wall and a small cavity. Thus the model helps us understand causes and consequences of pathology.
In the case of tricuspid valve hypoplasia, clinical decisions must be made early after birth. In a severe case of hypoplasia, the foramen ovale cannot be closed because of resulting unacceptably high pressure in the systemic veins caused by insuffiency of a severely hypotrophic right heart. Insertion of an aortic-pulmonary shunt is advised (12); however, the shunt has the large disadvantage that the right heart does not develop further. In the simulation with a moderate flow restriction shown in Fig. 8, this pressure is found to increase moderately. If no repair was made, foramen ovale flow was simulated to increase with age. Thus repair is predicted to be very useful, and one should not wait for more than a few weeks because of expected maladaptations. In less severe cases of flow restriction, foramen ovale flow was simulated to decline slowly to zero, meaning that no intervention is needed. The proposed decisions agree with clinical practice (12). Thus the CircAdapt model seems promising in predicting effects of interventions.
The sensitivity analysis (Table 4) shows adaptation and functional responses to various parameter changes. There are no extremely critical parameters, as shown by the moderate values in Table 4, all <250%. Chronic hypertension, as simulated by increase of mean aortic pressure (pref), resulted in concentric hypertrophy of the left ventricle, as expected. When the maximum level of volume load is increased chronically, as simulated by an increase of mean systemic flow during exercise, chamber wall volumes increase more than end-diastolic cavity volumes at rest, thus indicating eccentric hypertrophy. Left ventricular hypertrophy and stiffening of the passive matrix caused increase of pulmonary venous pressure in the simulation as well as experimentally (22).
Not all simulations appear realistic; e.g., at low atrial preload, atrial contractility appeared much higher in the simulation than in reality (Fig. 7). The aberrations were detected by inaccuracies in simulating pulmonary venous flow velocity. Because the CircAdapt model is largely based on undisputed physical principles, measured hemodynamics provide information about contractile properties of atrial tissue. Thus the model might be used to determine tissue properties in situ, which cannot be easily accessed in physiological experiments on isolated tissue.
The model has several limitations. We paid relatively little attention to adapting the physiological properties of the model to those of real physiology. All modules in the CircAdapt model may be improved or replaced, leaving the proposed structure unaffected. For instance, to simulate wave reflection and peripheral hemodynamics, modules simulating blood vessel segments with wave propagation must be developed. To simulate peripheral input impedances, windkessel models must be placed distally from the vessel network. For a more precise description, arteries and veins are likely to have different adaptation rules, and these rules may be age dependent (5, 18). The model of sarcomere mechanics has been derived from experiments on isolated cardiac muscle (14, 34). In such a preparation, the effect of cross-fiber stiffness is not included, whereas in the whole heart model, this stiffness is likely to be important, especially during diastole (39). Left-to-right ventricular interaction has been described by a coupling factor (19, 21, 27, 32). We modeled this interaction as an outer shell encapsulating an inner shell (Eq. 14) (1). In the fetal heart, the coupling was modeled to a symmetrical interaction. Modeling of ventricular interaction clearly requires further attention. For the atria, myocardial contractile characteristics are different (30), but more accurate data are not available. To obtain realistic behavior of mitral flow dynamics, atrial contractility was made more sensitive to stretch than ventricular myocardium (Eq. 20). In the valve module, the role of valve leaflet compliance, chordae tendinea, and papillary muscles has been neglected. Precise description of the hydrodynamic behavior of the valve requires numerical modeling of the flow field around the valves (25).
The modular setup makes the model flexible, as shown by the simulation of the fetal circulation. A further advantage of the modular setup is that a simple module can be replaced by a more complicated, e.g., a finite-element, model. Thus the CircAdapt model may provide realistic boundary conditions to the finite-element problem of simulating cardiac wall mechanics. The CircAdapt model is simple, needing relatively little computational effort. Because of flexibility, it has been chosen for simulation in MatLab, but in another programming language, e.g., C++, the model can be made much faster. Therefore, assessment of a therapy while you wait is a likely option.
In conclusion, the CircAdapt model realistically simulates total circulation dynamics and geometry. Combined application of a few, quite familiar rules, expressing adaptation of heart and blood vessels to mechanical load, suggests that the circulation forms as a result of self-structuring. Apparently, relatively little information, as condensed in a few adaptation rules, is sufficient to define the relatively complicated model of the circulatory system. Therefore, the number of unknown model parameters is much lower than in existing, more conventional models. Because most pathologies are isolated disorders, leaving the rest of the adaptation rules intact, the CircAdapt model is well suited by design to simulate individual pathologies with inherent adaptations. Also short- and long-term effects of pathology and therapeutic interventions can be predicted. The modular setup of the model enables easy simulation of aberrant circulatory pathology, e.g., occurring in congenital heart diseases.
Model of sarcomere mechanics.
According to Eq. 8, sarcomere (Cauchy fiber) stress (σf) and the derivatives of state variables were described as a function of time (t), sarcomere length (Ls), and the sarcomere-specific state variables. In the present simulations, we have used an empirical representation with a second-order differential equation, simulating experiments on isolated cardiac muscle (8) quite accurately; however, modifications were made to accommodate a lower heart rate, common to humans. Length of the contractile element (Lsc) and a time-variant contractility parameter (C) were state variables. Duration of contraction (tact) and contractility (Cs) were supplied separately.
Most parameter values are substituted in the equations below. For the normalized length of the series elastic element (LsNorm) and contractile element velocity (dLsc/dt) (A1) At the beginning of contraction, dC/dt is positive and determined by rise time (trise), Lsc, and Cs. Dependence on Lsc represents increase of contractility with sarcomere length. Decay occurs after a delay time proportional to tact and also decreases linearly with Lsc. A decay time constant (td) has been chosen so that decay time is about equal to rise time. Thus we have used with (A2) The function fRise(t/trise) represents a smooth positive unity amplitude pulse with duration proportional to trise. For the active stress component, it was used as follows (A3) Multiplication by (Lsc − 1.51 μm)0.5 was performed to simulate decrease of developed stress with short sarcomere length at the end of contraction. σf is the sum of an active and a passive component (A4) For stretch, the applied passive characteristic (Eq. A1.4) is numerically very close to the zero-pole model (14) but has the convenience of not having a mathematical singularity. For compression, passive myofiber stress is set to zero, thus simulating buckling of the myocardial wall. Coefficients σpas and σact were chosen, so that for Ls = 2.3 μm, passive stress (σf,pas) was 6.3 kPa and isometric active stress (σf,act) with Cs = 1 was 138 kPa, respectively (Table 1). For atrial myocardium, twitch duration is shorter (29), and related to that property, νmax is likely to be larger (∼1.5 times). Also, the range of sarcomere length has been chosen slightly larger to accommodate a large passive stretch at ventricular end systole and a normal range of sarcomere length for active contraction. Not as much quantitative information has been reported about atrial as about ventricular contractile properties.
This work was supported by Netherlands Heart Foundation Grant 2000T036 and University Hospital Maastricht Grant PF155. T. Delhaas is a Clinical Fellow of The Netherlands Heart Foundation (Dr. E. Dekker Fund).
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
- Copyright © 2005 by the American Physiological Society