Authors: Bettine G. van Willigen, M. Beatrijs van der Hout‐van der Jagt, Peter H. M. Bovendeerd, Wouter Huberts, Frans N. van de Vosse
Categories: Applied Research, Doppler ultrasound, fetal growth, mathematical modeling, pulsatility index, wave propagation
Source: International Journal for Numerical Methods in Biomedical Engineering
Doi: 10.1002/cnm.3877
Doppler ultrasound is a commonly used method to assess hemodynamics of the fetal cardiovascular system and to monitor the well‐being of the fetus. Indices based on the velocity profile are often used for diagnosis. However, precisely linking these indices to specific underlying physiology factors is challenging. Several influences, including wave reflections, fetal growth, vessel stiffness, and resistance distal to the vessel, contribute to these indices. Understanding these data is essential for making informed clinical decisions. Mathematical models can be used to investigate the relation between velocity profiles and physiological properties. This study presents a mathematical model designed to simulate velocity wave propagation throughout the fetal cardiovascular system, facilitating the assessment of factors influencing velocity‐based indices. The model combines a one‐fiber model of the heart with a 1D wave propagation model describing the larger vessels of the circulatory system and a lumped parameter model for the microcirculation. Fetal growth from 20 to 40 weeks of gestational age is incorporated by adjusting cardiac and circulatory parameter settings according to scaling laws. The model's results, including cardiac function, cardiac output distribution, and volume distribution, show a good agreement with literature studies for a growing healthy fetus from 20 to 40 weeks. In addition, Doppler indices are simulated in various vessels and agree with literature as well. In conclusion, this study introduces a novel closed‐loop 0D‐1D mathematical model that has been verified against literature studies. This model offers a valuable platform for analyzing factors influencing velocity‐based indices in the fetal cardiovascular system.
Keywords: Doppler ultrasound, fetal growth, mathematical modeling, pulsatility index, wave propagation
Doppler ultrasound is a commonly used method to assess hemodynamics of the fetal cardiovascular system and to monitor the well‐being of the fetus. Indices based on the velocity profile over time are often used for diagnosis, since they are expected to reflect the underlying physiology. For example, clinicians use the pulsatility index based on the velocity profiles over time of the umbilical and/or middle cerebral artery as an indirect indicator for placental resistance; an elevated placental resistance results in wave reflections causing an increase in pulsatility index [1], which is indicative for fetal growth restriction. Interpretation of these indices is difficult, since several factors influence these indices. With advancing gestational age, for instance, pulsatility index is affected by the relative increase of maximum versus minimum umbilical artery blood velocity [2]. The location of the measurement in the vessel impacts the pulsatility index, because velocity waveforms composed of forward and backward traveling velocity waves differ from shape and timing along the vessel's length [2, 3]. Additionally, placental resistance and vessel stiffness affect the forward and backward velocity waves, consequently influencing the pulsatility index as well [4]. These examples emphasize the significance of understanding the Doppler data to make well‐informed clinical decisions.
Mathematical models that incorporate fetal growth and the resulting changes in hemodynamics can provide better insights into (patho)physiology and can support clinical decision making by leveraging data interpretation [5, 6]. Such a model should include submodels of the fetal blood circulatory system and fetal growth.
The heart serves as the driving force for flow in the cardiovascular system. Most mathematical models that describe the fetal cardiac contraction use a time‐varying elastance model [7, 8, 9]. This model prescribes the contraction by varying the elastance of the cardiac chamber over time. Subsequently, fetal cardiac growth can be captured by scaling the minimal and maximal elastance [10]. This model has demonstrated to be able to simulate average healthy fetal cardiac growth realistically in terms of cardiac output, stroke volume, and ejection fraction. However, in pathophysiological situations, the elastance model is less suitable, because it fails to describe the cardiac adaptations due to fetal growth that lead to changes in hemodynamics and cardiac function. During fetal growth, the heart undergoes physiological adaptations, such as an increase of tissue stiffness and contractility, and increase of cavity and wall volume with advancing gestational age [11], while the time‐varying elastance model describes all those changes with one parameter, the elastance. Hence, simulating pathophysiology is more challenging. The one‐fiber model [12] describes cardiac function based on material properties and cardiac size and is, therefore, hypothesized to be able to capture fetal cardiac growth more realistically, especially when considering pathophysiological scenarios. Previous studies [13, 14] successfully scaled the one‐fiber model from adult to a full‐term fetus, but did not include fetal growth.
The current existing fetal mathematical models use a diode to characterize the cardiac valves [8, 14, 15], or incorporate variations on the Bernoulli equation [7, 13]. The model parameters in the latter approach are case‐specific and lack a unique model parameter set. Consequently, applying these valve models to different scenarios requires tuning and optimizing of the parameters. Hence, incorporating fetal cardiac growth is challenging. Therefore, in this study, the valves are described with diodes, allowing for fetal growth by scaling the valve resistance [10].
The circulatory system is often described using lumped parameter models [8, 13, 14, 16, 17]. These models are suitable to capture cardiac preload and afterload, but they fail to realistically describe propagation of velocity waves in the large vessels. However, velocity‐based Doppler indices depend critically on an accurate prediction of these waves. Therefore, 1D modeling is essential for better interpretation of diagnostic values, such as indices based on Doppler ultrasound measurements, during fetal development and different (patho)physiological scenarios [18].
Details of the microcirculation are often unknown, and the absence of significant wave propagation, due to high compliance and low flow, makes global simulations of flow, pressure, and volume sufficient for a general overview in the microcirculation. Additionally, these simulations establish boundary conditions for the wave propagation model describing larger vessels. Therefore, 0D modeling is allowed as this combines mathematical efficiency with sufficient clinical parameters.
Modeling of fetal growth is previously simulated with scaling laws in combination with lumped parameter models and time‐varying elastance models [10, 19]. However, fetal cardiovascular models that simulate fetal growth are limited and the combination with 1D modeling has not been made yet [6].
Therefore, this study aims to develop a mathematical model of the fetal blood circulation from 20 to 40 weeks of gestational age, focusing simulating velocity profiles for comparison with clinical Doppler ultrasound measurements. Verification of the model will be performed by comparing model outcomes with literature studies on healthy fetal growth to obtain a platform for analyzing factors that influence velocity‐based indices in the fetal cardiac system.
The overview of this study is as first, a mathematical description of the submodels (fetal blood circulation and fetal growth model) and the methodology of the Doppler indices simulations are given (Section 2). Subsequently, model outcomes, including cardiac function, cardiac distribution, and volume, are compared with literature data to verify the simulation of healthy fetal growth. After verifying this, Doppler indices are computed in various vessels and compared with literature data (Section 3). Finally, this study ends with a discussion about the model results, limitations, and an overall conclusion (Section 4).
The mathematical model introduced in this study includes two a fetal blood circulation model and a fetal growth model. The fetal blood circulation model integrates a one‐fiber model describing the cardiac function, a valve model replicating valvular motion, a 1D model capturing wave propagation through the larger vessels of the cardiovascular system, and a 0D model characterizing the microcirculation. The fetal growth model scales the parameters of the fetal blood circulation model from 20 to 40 weeks of gestational age. Figure 1 illustrates the interaction between these two models (elaborated in Section 2.2). The following sections give a mathematical description of both submodels, followed by the methodology that compares simulation results with Doppler ultrasound measurements.
FIGURE 1 A schematic representation of the interaction between two A fetal blood circulation and a fetal growth model. Fetal weight (w) and heart rate (fHR) undergo scaling based on gestational age (tGA in weeks) or are measured, when available, in the context of patient‐specific modeling. Subsequently, fetal weight and heart rate are used to scale full‐term fetal model parameters θw40fHR40 to new model parameters θwtGAfHRtGA that correspond to the specific gestational age, using the fetal growth model. These parameters are then applied to the fetal blood circulation model, resulting in pressure (p), flow (q), volume (V), and velocity (v) over time in each 0D and 1D element.
The contraction of each cardiac chamber is described using a one‐fiber model [12, 14]. This model relates the pressure inside a chamber (pcav) to the myocardial stress (σm), the cavity volume (Vcav), and the wall volume (Vwall) through the following
The myocardial stress σm includes passive stress σp (see Supporting Information: Equation A.3) and active stress σa. The active stress σa is expressed as the product of contractility c and three functions describing (1) the dependency on fiber length a1ls, (2) the time‐dependent activation function a2t, and (3) the fiber shortening dependency a3vs. Mathematically, the active stress is written as
The time‐dependent activation function a2t depends on the cardiac cycle time (T in seconds). The shape of the activation functions (Figure 2B) is derived from activation rise and decay times. These parameters are defined according to the equations provided by Mynard et al. [20] (see Supporting Information: Equations A.9 and A.10). Figure 2A illustrates the connection between the one‐fiber models (circles) describing the four heart chambers, the valves (Section 2.1.2), and the large vessels (Section 2.1.3). The pulmonary veins and vena cava connect to the left and right atrium, respectively, via a characteristic impedance (Zpv and Zvc, respectively). For a detailed explanation of the one‐fiber model and its model parameters describing a full‐term fetus, the reader is referred to Supporting Information: Appendix A.1.
FIGURE 2 (A) A schematic overview of the connection between the one‐fiber models (circles) describing the four heart chambers, valves (diodes and resistances), and large vessels. (B) The activation function (α2) of the ventricles and atria. AV: aortic valve, FO: foramen ovale, LA: left atrium, LV: left ventricle, MV: mitral valve, RA: right atrium, RV: right ventricle, TV: tricuspid valve, Zao: aortic characteristic impedance, Zpa: pulmonary artery characteristic impedance, Zpv: pulmonary venous characteristic impedance, Zvc: venous characteristic impedance.
The cardiac valves are described as ideal diodes. Hence, the valves have two states, completely open or closed, depending on the difference between the proximal pressure pp and distal pd pressure across the
where Rv,i represents the resistance of valve i, respectively, with i∈MVTVAVPVFO where MV the mitral valve, TV the tricuspid valve, AV the aortic valve, PV the pulmonary valve, and FO the foramen ovale. The resistance values representative for a full‐term fetus can be found in Supporting Information: Appendix A.2.
Pressure and flow wave propagation in the larger vessels is described by using 1D modeling (Figure 3). The 1D wave propagation of blood pressure (p) and flow (q) is governed by mass and momentum balance equations, assuming Newtonian incompressible fluid, a symmetric tube, and neglecting leakage (i.e., small side branches) and body forces (e.g., gravity forces). The approximate velocity profile [21] is used to express the wall shear stress and steady inertia forces in terms of p and q. Additionally, a constitutive relation for the arteries and veins is defined to describe the relationship between pressure p and cross sectional area using the area compliance [18]. The arteries are approximated as thin‐walled linear elastic tubes since their wall thickness‐to‐diameter ratio is smaller than 0.1. For the thin‐walled venous vessels, a power law that accounts for vessel collapse is used to describe their wall behavior [7]. For the umbilical arteries and vein, a thick‐walled linear elastic tube model is assumed due to their larger wall thickness‐to‐diameter ratio of 0.31 [22, 23]. For a detailed description of the mass and momentum balance, approximate velocity profile, constitutive relations, and estimation of model parameters for a full‐term fetus, the reader is referred to Supporting Information: Appendix A.3.
FIGURE 3 Schematic overview of the 1D elements describing the larger arteries (red) and veins (blue) of the fetus. This geometry is adopted from the study of Zhang, Haneishi, and Liu [9]. The labels refer to the vessel ID's. Matching the vessel ID's, the parameter values of the large vessels for a full‐term fetus can be found in Supporting Information: Table A.3. The connection with the microcirculation is represented in Figure 4.
The microcirculation of the fetal cardiovascular system is simplified using 0D elements consisting of resistances R, resistance impedances Z, and compliances C. Figure 4 shows a schematic overview of the microcirculation compartments. The labels correspond with the vessel ID's in Figure 3. The external pressure reflects the intra‐uterine pressure rather than the atmospheric pressure. It is identical for all vessels and set at pext=0 mmHg for all arterial Cart and venous Cven compliances. To avoid nonphysiological high frequency wave reflections between the 0D and 1D elements, characteristic impedances (Z) are added to connect the 0D and 1D
where Cj1D and Aref,j are the compliance and the cross‐sectional area at reference transmural pressure, respectively, of the connected 1D vessel j. For the methods to estimate the model parameters of the lumped parameter elements describing a full‐term fetus, the reader is referred to Supporting Information: Appendix A.4.
FIGURE 4 Schematic overview of the 0D elements describing microcirculation of the fetus consisting of compliances C, resistances R, and resistance impedances Z. The labels correspond with the vessel ID of the large vessels in Figure 3. Subscript abbreviations mean the art; arterial side, ven; venous side, PoV; portal vein, PS; portal sinus, ha; hepatic artery, liv; liver, hv; hepatic vein. This geometry is based on the study of Zhang, Haneishi, and Liu [9].
Most of the clinical data available for model verification relates to the healthy full‐term fetus. Consequently, the fetal blood circulation model was initially developed to describe a full‐term fetus. Typically, the parameter set θ, that depends on fetal weight w and heart rate fHR, to describe a full‐term fetus θw40fHR40 is either estimated from literature data or adopted from existing mathematical models. To extend the applicability of the model to describe fetal hemodynamics across gestational ages from 20 to 40 weeks, fetal weight and heart rate corresponding to each age are used to scale the parameter set from its original state at 40 weeks to a new parameter set θwtGAfHRtGA. Given that fetal weight and heart rate are routinely measured in clinical practice, this scaling approach facilitates the possibility of patient‐specific modeling. The output parameters of the fetal blood circulation model for each week between 20 and 40 weeks consists of pressure, flow, volume, and velocity within every 1D and 0D element included in the anatomical geometry (Figures 3 and 4). A schematic overview of the scaling methodology is illustrated in Figure 1.
In this study, global healthy fetal growth is assumed, and fetal weight and heart rate are determined based on their relation with gestational age tGA. The relation between gestational age and weight, based on the study of Kiserud et al. [24], can be described by the following
with w1, w2, w3, and w4 constants (Table 1) and t40=40 weeks. The relation between gestational age and heart rate, adopted from the study of Yigit et al. [10], can be expressed as
with f1 and f2 constants (Table 1). The cardiac cycle time (in seconds) is obtained by T=60fHR. Figure 5 gives a visual representation of Equations (5) and (6).
FIGURE 5 Visual representation of fetal weight and heart rate during fetal growth, ranging from 20 to 40 weeks gestational age (see Equations 5 and 6).
Model parameters for a full‐term fetus can be more accurately estimated, because of availability of literature data and models for that age group. Unfortunately, comparable information for lower fetal ages is not readily accessible. Therefore, the model parameter set of a full‐term fetus is scaled dependent on fetal weight and cardiac cycle time to estimate parameters for earlier ages. Sections 2.2.1, 2.2.4 elaborate submodel‐specific on which parameters are weight‐ or cardiac cycle‐dependent parameters. The weight‐dependent parameters Y are scaled with the following scaling
Note that Y40=Ywt40 is the corresponding parameter value for 40 weeks (see Supporting Information: Appendix A), and γ the scaling parameter. In addition, the model includes a set of cardiac cycle‐depending parameters X, given
where X40=XTt40 is the corresponding parameter value for 40 weeks (see Supporting Information: Appendix A). The following sections provide an overview of the scaling methodology for the model parameters for each fetal blood circulation submodel. Figure 6 represents the effect of the fetal growth model on the parameters within each submodel.
FIGURE 6 Visual representation of the normalized change in parameter values of each fetal blood circulation submodel, induced by the fetal growth model, from 20 to 40 weeks of gestational age (elaborated in Sections 2.2.1, 2.2.4). The dashed black line represents Y/Y40=1.
The weight‐dependent parameters of the one‐fiber model are presented in Table 2 with their scaling factor γ (Equation 7). The top left panel of Figure 6 visualizes the effect of the fetal growth model on the one‐fiber model. This paragraph elaborates on how the scaling factor is selected for the different parameters.
The association between the unloaded cardiac cavity volume Vcav,0 and fetal weight is assumed to be linear. The ventricular wall‐to‐unloaded cavity volume ratio rvol,v is estimated based on fetal cardiac dimensions of García‐Otero et al. [26] and is approximated to increase from 1.25 to 1.50 between 20 and 40 weeks, depending on fetal weight. In contrast, the atrial wall‐to‐unloaded cavity volume ratio rvol,a is assumed to decrease from 0.16 to 0.10 between 20 and 40 weeks, to maintain a relatively constant mean atrial pressure as demonstrated in Johnson et al. [27], while the cavity volume undergoes growth. To compensate for the pressure increase resulting from increased cavity volume, the arterial wall‐to‐unloaded cavity volume ratio is assumed to decrease. The specific decrease from 0.16 to 0.10 is tuned to ensure that the atrial pressure remains approximately constant.
In addition to these geometric adaptations of cavity and wall volumes based on fetal weight, the myocardium undergoes material changes as well. These material adaptations include alterations in contractility and fiber stiffness, leading to variations in contraction force. Contractility (c), an intrinsic characteristic of the myocardium, depends highly on the availability of and sensitivity to calcium and the number of sarcomeres, and it roughly doubles from 20 to 40 weeks [28, 29]. The scaling law for contractility c is defined to account for this doubling. Fiber stiffness, the ability of the myocardial fibers to stretch, depends on gestational age. Studies by van Oostrum et al. [30] and Elmstedt et al. [31] show reduced myocardial ventricular shortening with gestational age, indicating less elongation of fibers. In the model, this is simulated by decreasing the active stress scaling factor (σa,0) to mimic the weaker force of contraction due to the less elongation of the fibers with advancing gestational age. To simulate the increase in the ventricular compliance as demonstrated by Chang et al. [32], the passive stress scaling factor in the radial (σr,0) and fiber (σf,0) directions decrease with increasing fetal weight. The scaling factor is chosen such that a realistic E/A ratio of the atrioventricular valves is achieved. The sarcomeres are presumed to maintain uniform. Consequently, the curvature parameters (cf, cr, and ca) are assumed to remain constant throughout fetal growth, as well as the fiber lengths at zero passive stress (ls,0) and zero active stress (la,0).
The remaining parameters (vs,0, τr,0, τd,0, ar, ad, and t0) are cardiac cycle‐dependent parameters (see Equation 8).
The resistance of the atrioventricular valves and foramen ovale are weight‐dependent parameters, see Table 3, and the respective scaling factors are adopted from the study of Yigit et al. [10]. The top right panel of Figure 6 illustrates the effect of the fetal growth model on the resistances of the atrioventricular valves and the foramen ovale. For parameters describing a full‐term fetus, see Supporting Information: Appendix A.2. The resistances of the semilunar valves are determined based on the characteristic impedance of the connected 1D vessel (see Equation 4).
The parameters of the wave propagation model are weight‐dependent (Equation 7). The dynamic viscosity of blood η is associated with gestational age as described by Welch et al. [33]. This study translated this association to a dependency on fetal weight. The radius a0 and length l of the 1D vessels, as well as the wall thickness of the umbilical cord vessels hUC, undergo scaling through allometric scaling principles [34]. The scaling factor for the Young's modulus of the umbilical cord vessels and the stiffness parameter of the arteries Sa and Sv are chosen, such that the pulsatility indices are realistic for all gestational ages. The bottom left panel of Figure 6 shows the parameters of the 1D model, affected by the fetal growth model, from 20 to 40 weeks of gestation. For the reference values (40 weeks of gestational age), see Supporting Information: Appendix A.3.1. The scaling factors for the parameters of the wave propagation model are presented in Table 4.
The resistances, compliances, and unloaded volumes of the 0D model are scaled with fetal body weight (Equation 7). The scaling factor for the resistance is derived by estimating the resistance on 20 weeks based on a target flow q [35, 36] and pressure drop Δp [37] at 20 weeks R=Δpq. Subsequently, the scaling factor is defined, such that the scaling law agrees with the predefined values at 20 and 40 weeks. The arterial compliance is determined based on RC‐time τRC=RCart, which is scaled with fetal weight. The scaling factor of τRC is chosen, such that the maximal velocity of the arteries are realistic resulting in an increase from 0.16 at 20 weeks to 0.83 at 40 weeks. The venous compliance is 14.3 times larger than the arterial compliance, like the neonatal and adult circulation [7], and this compliance ratio is assumed to be constant during fetal growth.
The total unloaded volume is linearly scaled with fetal weight. The scaling law for the unloaded placental volume is chosen such that the placental blood volume fraction decreases from approximately 0.50 to 0.17 [38, 39] between 20 and 40 weeks of gestational age. The pulmonary unloaded blood volume fraction is assumed to remain 0.10. The rest of the unloaded volume is spread over the other 0D elements maintaining the volume ratio between the elements as for a full‐term fetus. The scaling factors for the parameters of the lumped parameter model are presented in Table 5. The bottom right panel of Figure 6 visualizes adaptation in resistance, volume and RC‐time of the microcirculation due to the fetal growth model. For full‐term fetus parameter values, see Supporting Information: Appendix A.4.
The solution of the fetal blood circulation computations involves an algorithm that consists of a cardiac cycle and hemodynamic convergence loop. In the first loop, pressure and flow are solved at each time step Δtc until the cardiac cycle time T is reached. Subsequently, the maximal difference between the pressure p and flow q of the current and previous cardiac cycles is determined using the normalized maximal root mean squared error ϵk:
Here, k represents the cardiac cycle number, t the time step within a cardiac cycle, and y can be replaced with either p or q. When ϵk is smaller than a predefined tolerance ϵtol (Table 6), hemodynamic steady state (i.e., a time periodic solution) is considered reached. Otherwise, the cardiac cycle loop restarts.
To initialize the simulation, mean circulatory mean filling pressure (pmcfp) is imposed on the entire fetal cardiovascular system and the total fetal blood volume (Vtot) is set to 110 mL/kg [40]. This initial pressure is based on the loaded volume (Vloaded=Vtot−Vunloaded) and the total compliance Ctot of the fetal cardiovascular
The loaded volume represents the volume that contributes to the pressure in the system above zero transmural pressure. The unloaded and initial volume in the 1D elements is based on the unloaded geometry parameters assuming a cylindrical tube. The unloaded volume of the 0D elements is obtained by literature and scaling. The initial volume in the 0D elements is computed from the unloaded volume of the element V0,e (see Supporting Information: Appendix A.4), the additional volume based on pmcfp and the corresponding compliance Ce, according
To couple the 0D and 1D elements, the method described by the study of Kroon et al. [41] is employed. This method reformulates the wave propagation and lumped parameter model equations into a unified
which allows the computation of the pressure p¯t+Δt and flow q¯t+Δt vectors for the next time step (length N) based on the pressure and flow of previous time steps using a second‐order backward scheme. The information from previous time steps is collected in the right side vector f¯ of length N. The matrix K of size N×N contains the characteristics of the elements. To solve the pressure for the next time step p¯t+Δt, the inflow and outflow of an element are assumed to be in opposite directions resulting in a q¯t+Δt=0 for all nonboundary nodes. The boundaries nodes are either prescribed with an external pressure pex=0 or connected to a one‐fiber model. Hence, each node has pressure or flow information, which allows us to solve pressure p¯t+Δt. Subsequently, flow for the next time step q¯t+Δt can be computed using p¯t+Δt.
The cardiac function model operates as a decoupled submodel, meaning that the pressure within the connected 1D vessels is used to compute the in‐ and outflow of the cardiac chambers. These flows are then incorporated into the vector f as external flows. Simulation parameters can be found in Table 6.
The source code is written in Python and is available on request on Zenodo ^1^ .
Typically, the maximal envelope of the Doppler velocity profiles is used to compute Doppler indices, that is, the maximal velocity in the Doppler sample volume. Therefore, to allow model verification with clinical data, simulation results on the maximal velocity in a vessel are required. To simulate the maximal velocity, a velocity profile has to be assumed. Typically, a factor h is multiplied with the mean velocity (v^zz,t) to obtain the maximal velocity maxvzz,t=hv^zz,t with h=2 for a fully developed velocity profile and h=1 for a flat velocity profile. However, this factor implies that the profile is constant for the entire vessel and also stable over time. In this study, the approximate velocity profile is used, which takes time and location into account [21]. This velocity profile divides the vessel into a central core where viscous forces are negligible and a viscous layer where friction and pressure gradient are the main forces acting on the fluid. The velocity profile vzr,z,t on location z along the radius a (perpendicular to the flow direction [z]) and time t point is then given
with the mean cross‐sectional velocity v^zz,t on location z, η the dynamic blood viscosity, and ∂p∂z the pressure gradient over the segment. The functions ϕ1 and ϕ2 are defined as
with ζ^=maxart2ζc and ζc=max0,1−2α with α the Womersly number and
To convert the mean cross‐sectional velocity v^zz,t to maximal flow velocity, a conversion factor for every vessel is computed as
Note that the mean cross‐sectional velocity v^zz,t is obtained by the model based on the simulated flow qzz,t and cross‐sectional area Az,t v^zz,t=qzz,tAz,t.
To verify the ability of the mathematical model to replicate velocity‐based Doppler ultrasound indices during fetal development, a systematical analysis of the model simulations is performed. Healthy fetal growth is simulated in 4‐week intervals from 20 to 40 weeks of gestation. Subsequently, an overall analysis of fetal growth is conducted, followed by an examination of velocity profiles and computed Doppler indices. Next sections delve into the specifics of the analysis methods.
Doppler ultrasound indices depend on various factors, including fetal growth. Therefore, before verifying Doppler indices computations, it is essential to confirm the realism of healthy fetal growth simulations. Given the difference between the fetal and adult cardiovascular system, specific cardiac function parameter values are characteristic of fetal circulation. To ensure realistic healthy fetal growth, the simulated cardiac function parameters are verified against literature data. These parameters include left (LCO), right (RCO), and combined (CCO) cardiac output, left‐to‐right cardiac output ratio (RCO/LCO), right (REF) and left (LEF) stroke volume, right (REF) and left (LEF) ejection fraction, and mitral and tricuspid valve E/A ratio. Additionally, an analysis of the distribution of cardiac output and blood volume over fetal growth is conducted, with comparisons to available data. For the results, the reader is referred to Section 3.1.
Following the assessment of healthy fetal growth, velocity profiles are evaluated by comparing the peak systolic velocity vPSV, the pulsatility index IPI (the ratio between the difference between the systolic and diastolic velocity and the mean velocity IPI=vPSV−vPDVv^zz,t), and the pulsatility index for veins IPIV (the ratio between the difference between the systolic and atrial velocity and mean velocity IPIV=vPDV−vAv^zz,t). First, an analysis of velocity wave profiles in the ascending aortic root and umbilical artery, from core to vessel wall, is conducted to demonstrate the factors influencing the velocity profile. This analysis also serves to confirm the model's capability to simulate the variation of the pulsatility index IPI along the length of the umbilical artery. Subsequently, vPSV, IPI, or IPIV in various vessels are computed for different gestational ages and compared with literature. For the results, the reader is referred to Section 3.2.
Figure 7 illustrates the cardiac function of the simulated growing healthy fetus from 20 to 40 weeks (corresponding to fetal weight increase from 0.31 to 3.66 kg [24]) compared with regression lines or mean values with standard deviation based on literature data [42, 43, 44, 45, 46, 47, 48]. The simulated left (LCO) and right cardiac output (RCO) increase from 77 and 91 mL/min at 20 weeks to 607 and 784 mL/min at 40 weeks, respectively, resulting in an combined cardiac output (CCO) increase from 168 to 1391 mL/min. Throughout gestational age, the simulated right‐to‐left cardiac output ratio (RCO/LCO) remains approximately 1.25, indicating a right dominant heart. Additionally, the simulated left (LSV) and right (RSV) stroke volume increase from 0.51 and 0.61 mL at 20 weeks to 4.3 and 5.6 mL at 40 weeks, respectively. The simulated left (LEF) and right (REF) ejection fraction decrease from 0.66 and 0.63 at 20 weeks to 0.59 and 0.57 at 40 weeks, respectively. The simulated E/A ratio of the mitral and tricuspid valve increase from 0.57 and 0.64 at 20 weeks to 0.72 and 0.76 at 40 weeks. These simulated results show the same trends and realistic values for cardiac function as the data from literature (Figure 7).
FIGURE 7 Overview of simulated (green diamonds) cardiac function output parameters illustrating global healthy fetal cardiac growth compared with regression lines (black dashed lines) or mean values with standard deviation (black marker and hat) based on literature data. *Median and 10% and 90% confidence interval.
Figure 8 presents the mean flow to the organ regions illustrated in Figure 4 as well as within the fetal shunts, both compared with mean flow measurements with standard deviations derived from magnetic resonance imaging techniques [36] or Doppler ultrasound measurements [35, 39]. The simulated mean flow values with advancing gestational age fall within the standard deviation reported in the literature data, with the exception of a slightly lower flow through the ductus arteriosus at 40 weeks. This discrepancy leads to a reduced flow to the upper body, falling below the standard deviation. The cardiac output distribution increases to the liver from 20% at 20 weeks to 28% at 40 weeks. The study of Rudolph [49] reports a cardiac distribution of 25.2% to the liver in late‐gestation lambs. The cardiac output distribution to the lower body decreases from 17% at 20 weeks to 11% at 40 weeks. The percentage of shunting from the umbilical vein to the ductus venosus (DV) remains about 29%. The study of Kiserud and Acharya [39] reports a shunting of 28% at 20 weeks and approximately 20% from 24 weeks onwards.
FIGURE 8 Flow to different organ regions and in the fetal shunts compared with literature [35, 36, 39] (mean ± standard deviation). DA: ductus arteriosus, DV: ductus venosus, and FO: foramen ovale.
Figure 9 shows the simulated volume distribution to the placental and fetal circulation. The figure shows an increase of blood volume within the fetal circulation relatively to the total blood volume from 44% at 20 weeks to 77% at 40 weeks of gestation. Consequently, the volume within the placental circulation decreases from 56% to 23%. These findings follow the trend described in the research of Kiserud and Acharya [39]. The blood volume in the pulmonary circulation remains relatively stable at 11%.
FIGURE 9 Volume distribution to placental and fetal circulation from 20 to 40 weeks compared with literature data [39].
Figure 10 illustrates the influence of specific time point and measurement location in the cardiac cycle on the velocity profile. The upper figures show velocity and pressure gradient information of the ascending aorta located just after the aortic valve (Womersley number α=8.74). In the left panel, the velocity profile is shown across the cross‐section from the core (radius = 0 mm) to the maximal radius, during various time points of the cardiac cycle (colors). The maximal cross‐sectional velocity varies over time. In addition, the shape of the velocity profile within the viscous boundary layer varies over time with the pressure gradient dpdz (right panel). A negative pressure gradient results in a higher velocity than within the core, while a positive pressure gradient leads to a lower or even negative velocity value. In addition, the right panel shows the maximal cross‐sectional velocity throughout a cardiac cycle (representing the maximal envelope of the Doppler ultrasound velocity profile) with the time points colored, mirroring the left panel.
FIGURE 10 The upper figures represent the velocity profiles and pressure gradient of the ascending aorta close to the aortic valve. Lower figures show the same information for the umbilical artery at fetal end (straight line) and placental end (dashed line). The left panel demonstrates the velocity profile of the corresponding vessel across the cross section of the vessel, from the core (radius = 0 mm) to the maximal radius, during different time points (colors). The right panel presents the maximal velocity throughout a cardiac cycle at fetal and placental end with the time points highlighted, mirroring the left panel. In addition, the right panel contains the pressure gradient of the vessel's length dpdz.
The lower figures show the same information for the umbilical artery at two fetal end (straight line) and placental end (dashed line), demonstrating that the velocity profile depends on z location as well, besides time point and radial location. As the Womersley number of the umbilical artery is lower (α=4.42), the viscous boundary layer is longer with respect to the length of the vessel [21], demonstrating the influence of the size of the radius on the velocity profile. The right panel shows that the pulsatility index IPI varies with the location within the vessel with a IPI of 1.05 at the fetal end and 0.87 at the placental end. These findings align with values reported by Acharya et al. [2], demonstrating an agreement within 90% confidence interval of [0.58, 1.11] and [0.57, 0.99], respectively.
Figure 11 displays the maximal cross‐sectional velocity profiles at the origin of the ascending aorta, main pulmonary artery, left renal artery, and splenic artery, as well as at the center of the left umbilical artery (upper graphs). The lower graphs show the peak systolic velocity vPSV or pulsatility index IPI in comparison with literature data (mean ± CI [10%, 90%]). The vPSV in the aortic and main pulmonary artery increases from 0.69 and 0.64 m/s at 20 weeks to 0.86 and 0.90 m/s at 40 weeks, respectively, agreeing with literature data [48, 51]. Both velocity profiles present a clean peak without any notches or backflow.
FIGURE 11 The velocity profiles of the maximal cross‐sectional velocity of the ascending aorta (vessel ID = 1), main pulmonary artery (vessel ID = 46), left renal artery (vessel ID = 72), splenic artery (vessel ID = 24), and left umbilical artery (vessel ID = 45), as calculated with the model (upper graphs). Lower graphs show the model results for peak systolic velocity (vPSV) or pulsatility index (IPI) compared with literature data [2, 48, 50, 51, 52] (mean ± CI [10%, 90%]).
The simulated renal and splenic artery vPSV increase from 0.29 and 0.20 m/s at 20 weeks to 0.53 and 0.53 at 40 weeks, respectively. The vPSV of both vessels agree with Contag et al. [50] and Tongsong et al. [52]. The model simulation shows backflow during ventricular diastole for both renal and splenic artery.
The umbilical artery velocity profile resembles a skewed sinus wave profile, with minimal and maximal velocity value of 0.07 and 0.33 at 20 weeks and 0.19 and 0.44 m/s at 40 weeks, respectively. The study of [2] indeed indicates a fairly constant vPSV with advancing gestational age and an increase end‐diastolic velocity intra‐abdominally. The PI falls within the 90% confidence interval.
The venous velocity profiles of the superior vena cava outlet, inferior vena cava outlet, pulmonary vein outlet, and DV inlet adopt an M‐shaped pattern starting with a valley caused by atrial contraction vA followed by a systolic vPSV and diastolic vPDV peak (Figure 12). The pulsatility index for veins (IPIV) for these vessels is compared with literature data [25, 53, 54, 55] (mean ± CI [10%, 90%]) can be seen in the lower graphs. The M‐shape pattern agrees with the literature data and the IPIV of the vessels fall within the 90% confidence interval.
FIGURE 12 The velocity profiles of the maximal velocity of the superior vena cava (vessel ID = 50), inferior vena cava (vessel ID = 70), pulmonary vein (vessel ID = 82), and ductus venosus (vessel ID = 69) calculated with the model (upper graphs). Lower graphs show model results for the pulsatility index for veins (IPIV) for these vessels compared with literature data [25, 53, 54, 55] (mean ± CI [10%, 90%]).
To enhance interpretation of medical data, this study introduced a novel closed‐loop 0D‐1D mathematical model describing global healthy fetal growth with the emphasis on Doppler ultrasound measurements. These following paragraphs discuss major findings and limitations.
The available literature data that describe cardiac function display a wide range of values both within and between the datasets (Figure 7). Despite this variability, there is a general consensus regarding the overall trend observed in the datasets, with the exception of the right‐to‐left cardiac output ratio. The simulation results align with the clinically observed trends, and the right‐to‐left cardiac output ratio from the model falls within the range of the literature values. Additionally, the model simulates realistic flow to different organ regions (Figure 8) and distribution of volume (Figure 9). The model percentage of shunting from the umbilical vein to the DV remains about 29% throughout gestation from 20 to 40 weeks, while the study of Kiserud and Acharya [39] reports a shunting of 28% at 20 weeks and approximately 20% from 24 weeks onwards. This study shows a wide variation of results as well. Hence, the model results are considered acceptable. In other words, these results show that the model is capable of simulating global fetal growth during the second half of pregnancy. However, as the variation between and within the datasets is wide, patient‐specific modeling requires an in‐depth analysis of the importance of the model parameters.
During Doppler ultrasound measurements, a specific location within the vessel (Doppler sample volume) is selected to measure the velocity. Typically, the aim is to target the center of the vessel to capture the cross‐sectional maximal velocity, and two common assumptions for the velocity profile are a flat (maxvz=v^z) or a Poiseuille (maxvz=12v^z) velocity profile. Although the choice of velocity profile minimally impacts the analysis of pulsatility indices, it does influence the absolute velocity values, consequently affecting flow computation. In this model, flow is computed by assuming a velocity profile based on the assumption of a Stokes layer, which then needs to be translated back to cross‐sectional mean velocity. Dividing the flow by the vessel area yields the cross‐sectional mean velocity. Hence, for a meaningful comparison with Doppler measurements, the accurate definition of the velocity profile is crucial to translate the cross‐sectional mean into maximal velocity.
In this study, the approximate velocity profile [21] is used to define a velocity profile across the cross‐section of the vessel. This profile lies between flat and Poiseuille depending on the radius (i.e., reflected by Womersley number α) of the vessel. Figures 11 and 12 demonstrate good agreement with measured data, thus validating the acceptability of the chosen velocity profile. To highlight the importance of understanding the data before drawing diagnostic conclusions, Figure 13 displays the cardiac cycle mean of the ratio between cross‐sectional maximal and mean velocity fc against the radius of the vessel considered. As shown in the figure, Poiseuille profile is apparent in the smallest vessels, while flat profile is not reached in the largest vessels. In addition, Figure 10 demonstrates the difference in the shape of the viscous boundary layer and the maximal cross‐sectional velocity, both influenced by the pressure gradient and the Womersley number α of the vessel. These figures emphasize that the timing and location of Doppler measurements impact the measurement outcome. Such insights, as demonstrated in Figure 10, can be used to translate the clinical measured Doppler velocity profile to blood flow.
FIGURE 13 The cardiac cycle median of the ratio between the maximal and mean cross‐sectional velocity fc against the vessel radius.
There is limited information regarding the material properties of the fetal cardiac chambers with advancing gestational age. Therefore, in this study, it is presumed that the distinction in contraction force between atria and ventricles arises from the difference between cavity and wall volume. Consequently, it is assumed that the material properties of the ventricles and atria are equal. While the material composition may influence the contraction force, the thickness of the cardiac wall can compensate for varying stiffness levels. This rationale justifies the selection of these model parameter values and it reduces the number of unknown parameters. Similarly, due to scarce available information on the passive stress, it is assumed that the scaling factors for the passive stress in both the fiber and radial directions are equal. These model parameters were chosen to achieve realistic intracardiac pressures, thereby addressing potential differences in fiber and radial direction and minimizing the number of unknown parameters. Following the data from Johnson et al. [27], the atrial wall‐to‐cavity volume ratio is assumed to decrease from 0.16 to 0.10 between 20 and 40 weeks of gestational age, in order to maintain a relatively constant mean atrial pressure. However, they measured the right and left atrial pressures of only 6 and 13 fetuses, respectively, within the gestational age range of 19–30 weeks. As a result, the trend of fetal atrial pressures between 20 and 40 weeks of gestational age remains uncertain. Additional clinical measurements are required to investigate the impact of atrial pressure on flow and pressure dynamics during fetal growth.
The cardiac valves are described as diodes (Section 2.1.2), which conveniently do not need adaptation during fetal growth and result in realistic mitral and tricuspid valve E/A ratios (see Figure 7). However, to avoid oscillations, a resistance equal to the impedance of the connected 1D vessel is required. These resistances result in a substantial pressure drop with open valves, reaching approximately 9 mmHg at 20 weeks and 14 mmHg at 40 weeks. In adults, a pressure drop over the aortic valve of less than 5 mmHg is considered normal and less than 25 mmHg is categorized as mild stenotic [56, 57]. Since the pressure drop over fetal valves is unknown, it is expected that similar values apply for to fetal valves. However, the higher than expected pressure drop may lead to supraphysiological intracardiac pressures when realistic mean arterial pressure is obtained. The study of Leyh et al. [58] shows a three‐phase valve motion for the adult valve. Hence, the two‐phase representation (completely closed and open) does not accurately capture the opening and closing motion of the leaflets [58]. Consequently, these valves lead to nonphysiological flow patterns without notches or backflow (see Figure 11). Other fetal cardiovascular modeling studies use the valve model proposed by Mynard et al. [20], which describes valve behavior using Bernoulli equation combined with an opening and closing curve. While this model results in a more realistic valve behavior, it is tuned to a particular case and lacks a unique model parameter set. Therefore, applying this valve model to different scenarios requires parameter tuning and optimization. Hence, valvular growth cannot be incorporated by simply applying scaling laws.
Thus, although the instant opening and closing of the diodes may introduce nonphysiological wave reflections, affecting pressure and flow solutions, the currently available more sophisticated valve models require extensive tuning, which is not feasible when accounting for the dynamic changes during fetal growth. Additionally, it is likely that the diode models might have less pronounced effects in the fetal circulation due to presence of shunts. Improving valve modeling remains beneficial for capturing the complex dynamics of the fetal cardiovascular system and requires further research to minimize pressure drop and achieve a more realistic wave pattern throughout the cardiac system.
Venous valves play a role in ensuring unidirectional blood flow, particularly in response to changes in fetal position and pressure variations. However, including venous valves would increase the complexity of the model, especially given the limited detailed anatomical and physiological information available. Without precise data, addition of venous valves introduces more uncertainties. Fetal circulation is predominantly influenced by the heart valves and shunts, and changes in fetal position is not (yet) included in this model. Therefore, the main flow patterns are represented, and exclusion of venous valves is considered sufficient at this stage of the model development.
The dimensions of the 1D vessels in the model are based on Doppler ultrasound measurements obtained from literature. These measurements are used to determine the radius and area at unloaded pressure. However, it should be noted that during the Doppler ultrasound measurements, internal blood pressure is applied to the vessel. To address this issue, it is necessary to choose the appropriate reference pressure in the constitutive relation described in Supporting Information: Equation A.26. The reference pressure should correspond to the timing of the Doppler ultrasound measurements. Unfortunately, the specific timing is unknown in this study. In this model, diastolic pressure is selected as the reference pressure, considering that the diastolic phase is typically longer than the systolic phase. However, it is important to acknowledge that using the diastolic pressure as the reference may lead to an overestimation of the vessel dimensions and consequently an overestimation of vessel expansion.
The RC‐time constant and the ratio between venous and arterial compliance for a full‐term fetus are based on neonatal data and assumed equal for all regions. The unloaded volume is chosen, such that pmcfp=8 mmHg at 40 weeks, as this results in a mean atrial pressure of approximately 3.5 mmHg [27]. The mean circulatory filling pressure pmcfp is the resting pressure after the heart stops and blood is redistributed over the cardiovascular system. After applying scaling laws, pmcfp increased from 6 mmHg at 20 weeks to 8 mmHg at 40 weeks. For adults, pmcfp is around 7 mmHg [59].
The model simulates generic healthy fetal cardiac growth by scaling the input parameters of fetal weight and heart rate according to predefined growth in gestational age. These parameters are routinely measured in the clinic, which holds the potential for the model to support (patient‐specific) clinical decision‐making in the future. For the cardiac function model, the scaling laws are mostly defined by physiological reasoning based on information from literature, while the scaling laws for the 1D and 0D elements are based on previous modeling studies, on allometric scaling principles, or tuned to obtain physiological outcome. The results show realistic fetal cardiac growth, but the uniqueness of the parameters has not been established, nor their significance. To uncover the true underlying phenomena of growth, future research should investigate those two aspects.
Besides enhancing our knowledge about the fetal circulation, this novel mathematical model could serve as the foundation for a patient‐specific digital twin model of the fetal cardiovascular system, potentially enhancing clinical decision‐making. To achieve this, it is crucial to identify the most influential model parameters on simulation outcomes. These parameters should ideally be coupled in real‐time to measurements from the fetus, such as Doppler ultrasound measurements or fetal electrocardiograms. In cases where direct measurements of these parameters are unavailable or impractical, they can be estimated from comparison of model predictions with actual fetal outcomes. For instance, while a fetal electrocardiogram provides direct information about heart rate, placental resistance must be estimated from other measurements. Ultimately, digital twin technology offers additional insights that are otherwise unmeasurable, thereby supporting treatment planning through more informed diagnosis.
In conclusion, this study introduced a novel closed‐loop 0D‐1D mathematical model describing global healthy fetal growth to enhance the comprehension and interpretation of medical data with the emphasis on Doppler ultrasound measurements. Additionally, this model can potentially provide better insight into (patho)physiology and can support clinical decision making by leveraging data interpretation. It can also lead to serve as foundation for a digital twin and clinical decision support tools [6, 60]. This study introduces a parameter set and scaling laws that lead to realistic fetal growth and Doppler indices. However, the uniqueness or sensitivity of these parameters is not investigated. To uncover the underlying phenomena of growth, future research should investigate those two aspects.
B.G.W. wrote the main manuscript text based on a discussion session with all other authors. All authors contributed to the article and approved the submitted version.
The authors have nothing to report.
We declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
The data that support the findings of this study are available from the corresponding author upon reasonable request.