Closed-Loop Multiscale Computational Model of Human Blood Circulation. Applications to Ballistocardiography

Cardiac mechanical activity leads to periodic changes in the distribution of blood throughout the body, which causes micro-oscillations of the body’s center of mass and can be measured by ballistocardiography (BCG). However, many of the BCG findings are based on parameters whose origins are poorly understood. Here, we generate simulated multidimensional BCG signals based on a more exhaustive and accurate computational model of blood circulation than previous attempts. This model consists in a closed loop 0D-1D multiscale representation of the human blood circulation. The 0D elements include the cardiac chambers, cardiac valves, arterioles, capillaries, venules, and veins, while the 1D elements include 55 systemic and 57 pulmonary arteries. The simulated multidimensional BCG signal is computed based on the distribution of blood in the different compartments and their anatomical position given by whole-body magnetic resonance angiography on a healthy young subject. We use this model to analyze the elements affecting the BCG signal on its different axes, allowing a better interpretation of clinical records. We also evaluate the impact of filtering and healthy aging on the BCG signal. The results offer a better view of the physiological meaning of BCG, as compared to previous models considering mainly the contribution of the aorta and focusing on longitudinal acceleration BCG. The shape of experimental BCG signals can be reproduced, and their amplitudes are in the range of experimental records. The contributions of the cardiac chambers and the pulmonary circulation are non-negligible, especially on the lateral and transversal components of the velocity BCG signal. The shapes and amplitudes of the BCG waveforms are changing with age, and we propose a scaling law to estimate the pulse wave velocity based on the time intervals between the peaks of the acceleration BCG signal. We also suggest new formulas to estimate the stroke volume and its changes based on the BCG signal expressed in terms of acceleration and kinetic energy.


Normalized time-varying function in the model of cardiac elastance
In Equation 2, defining the time-varying elastance of a cardiac chamber, we use a normalized time varying function ( ) defined as in (Liang et al., 2009). Considering that the left and right ventricles (CC = LV and RV, respectively) contract simultaneously at the very beginning of the cardiac cycle, the latter function can be written: with the duration of the cardiac cycle, the duration of the ventricular contraction, and the duration of the ventricular relaxation. A similar model is chosen for the left and right atria (CC = LA and RA, respectively), which are also assumed to contract simultaneously. However, the time dependence of the corresponding normalized elasticity functions ( ) = ( ) is changed so that the atrial contraction starts at 0.8 (Gallo et al., 2020). We write the duration of the atrial contraction and the duration of the atrial relaxation. This gives the following equation for the atria: The numerical values used to characterize the time-varying elastance in each cardiac chamber are taken from (Liang et al., 2009) and given in Table 1 (core of the manuscript) for the case of normal contractility and = 0.86 (corresponding to a heart rate of 70 bpm).

Empirical laws used for the arteries of the healthy young case
In the systemic circulation, an empirical law is used to compute the pulse wave velocity 0 ( ) based on the radius ℛ( ) extracted from the segmentation of the systemic arterial tree: with the radius ℛ( ) measured in meters and the coefficients 1 = 4.0 × 10 5 kg. s −2 . m −1 , 2 = −500 m −1 , and 3 = 40.0 × 10 3 kg. s −2 . m −1 adapted from (Mynard and Smolich, 2015) to achieve normal wave velocities in both the large and smaller vessels (5.08 to 13.58 m. s −1 ).

Correction of the friction force at the wall in the case of non-established flow
In the case of a Poiseuille flow (i.e., a parabolic velocity profile) of an incompressible viscous fluid developed through the entire length of a straight cylindrical artery , the resistance to the flow is minimal and given by (Pedley et al., 1970a): with and ℛ the length and radius of the artery , respectively, and the dynamic viscosity of blood.
However, it has been shown that the hypothesis of a fully developed Poiseuille flow is not always valid, especially in the large proximal systemic arteries (Seed and Wood, 1971;Bogren and Buonocore, 1999;Tortoli et al., 2002). Hence, it has been suggested to correct the Poiseuille resistance to the flow by multiplying it by a factor (Pedley et al., 1970a(Pedley et al., , 1970b) that can be expressed as the ratio of the wall shear stress in the artery to the wall shear stress that would take place if the Poiseuille flow was established (see supplementary material in (Buess et al., 2020)): with ̅ the time average of the mean cross-sectional velocity in artery , and the volumetric mass density of blood.
If the evaluation of gives a result smaller than 1, it is then set to 1, in agreement with the fact that the Poiseuille resistance corresponds to the minimal resistance to blood flow. This means that we can introduce a critical length of the artery , , needed to establish a Poiseuille flow: This definition indeed leads to: giving > 1 when < , .
We can also express the non-establishment of the flow with the parameters and = +2 +1 defining the velocity profile in the artery , as is the ratio of the actual friction force at the wall per unit length ( , see Equation 11) to its value if the velocity profile was parabolic (i.e., if = 2): leading to: Based on this, an iterative approach is chosen to solve the equations defining blood flow in the arterial tree. First, we assume a Poiseuille flow (i.e., = 1) in all the arteries, which allows a first estimation of blood flow in each artery and the computation of the associated correction factors (or, alternatively, Coriolis coefficients). Then, the system of equations is solved a second time, considering the non-establishment of blood flow in some arteries.

Comparison with experimental results
The comparison of an artificial 3D BCG signal produced by a numerical model with an experimental 3D BCG signal is a relatively difficult task. Indeed, the signals along the and axes are very sensitive to the experimental conditions and, since their amplitudes are usually lower than the one of the y axis, they are easily contaminated by other phenomena (lower signal-to-noise ratios). Depending on the way BCG is recorded, these phenomena can include the precordial movements (contamination by the seismocardiographic signals, especially on the axis), the movements induced by the rotation of the body (especially on the axis (Baan et al., 1966)), the movements induced by local flow of blood below the BCG sensor (usually on the axis), or the movements that happen on the y axis but that are recorded on the or axes (depending on the exact position of the subject and of the BCG sensor).
Supplementary Figure 1 displays the 3D acceleration BCG signal generated by the numerical model in the young and healthy reference case and two types of experimental records: using a portable BCG device ( , , and ) attached to the lumbar region of a subject and using an ultra-low frequency BCG table ( and ) (van der Linde and Knoop, 1966). Except for the absence of a first upward H wave on the artificial signal, the agreement between these three BCG signals on the axis is very good, with a similar amplitude and position of the different waves. We also find a relatively similar amplitude on the axis of the three signals, but it is much more subject to noise, even though some similar patterns are found between the numerical model and the ultra-low frequency table. Finally, on the axis, while there is no ultra-low frequency BCG table data available for comparison, the portable BCG signal seems to be highly impacted by a phenomenon that is not considered in the numerical model. We assume that the large waves observed on the axis of the portable BCG signal correspond either to a contamination by the seismocardiographic signal or, most probably, by the propagation of a pulse wave through a vessel located just below the BCG sensor. It is thus very likely that this discrepancy between the numerical model and the portable BCG records are the consequence of some particularities of the portable BCG technique rather than a misestimation of the movement of the center of mass induced by cardiovascular activity.
Supplementary Figure 1. Comparison of the 3D acceleration BCG signal given by the numerical model with typical experimental signals given by a portable BCG device ( , , and ) attached to the lumbar region of a subject and by an ultra-low frequency BCG table ( and ). Data for the ultra-low frequency BCG table are taken from (van der Linde and Knoop, 1966).

6
Complementary results related to pulse wave velocity

Determination of the scaling law
Based on Equations 15 and A4, and the fact that the I, J, and K peaks are mainly caused by the propagation of blood flow along the aorta, we can write that in the aorta: with 0,Ao the average pulse wave velocity along the aorta and 0,Ao the average cross-sectional area at pressure ref along the aorta.
Between two peaks of the longitudinal signal, separated by a time interval ∆ , we assume that the slope of the pressure wave is proportional to 1 , and thus to 1 √ . We also assume that the difference in the amplitude of blood flow between two sections of the aorta is proportional to 1 .
Therefore, we can write: with ∆ measured in s 1/4 .cm -1 and being equivalent to a pulse transit time after a transformation operated on the ∆ time interval between two peaks on the BCG signal, while is a proportionality factor measured in s 3/4 .

Effect of the scaling law
The dataset described in the main text is used to assess the performance of the suggested scaling law (Equation 33) to evaluate the inverse of the mean pulse wave velocity in the aorta. Supplementary Figure 2 presents the results of the correlation analysis based on the time intervals directly measured on the longitudinal axis of after low-pass filtering with a cutoff frequency at 25 Hz. The results based on the same time intervals, but transformed with the suggested scaling law, are already presented in the main text (Figure 9).

Supplementary Figure 2.
Correlation of the inverse of pulse wave velocity in the aorta with a selection of time intervals: a. RJ; b. RK; c. IJ; and d. IK. All the results are obtained after low-pass filtering of the acceleration signal at 25 Hz. In each case, the global linear regression line is indicated together with its equation and the regression coefficient R 2 . Colors indicate different simulated ages, while symbols indicate different stressed volumes. Several points can have the same color and symbol, when they correspond to different heart rates. The healthy young reference case is indicated by a black diamond. Figure 2 display very bad correlation coefficients with the inverse of the pulse wave velocity in the aorta. These results differ a lot from those presented in Figure 9, for which the scaling law is used (e.g., R 2 = 0.02 for IK vs. R 2 = 0.98 for ̃) . The only exception is the RJ time interval, since its correlation coefficient is already very high without the scaling law (R 2 = 0.99). Similar results are also observed for the RI time interval (not presented).

Most of the time intervals presented in Supplementary
These results overall support the interest of using the scaling law to evaluate pulse wave velocity in the aorta.

Effect of filtering
It is reasonable to assume that different body types will have different dampening effects on the BCG signal. The sensitivity of a metric to these dampening effects can be studied by looking at different filters applied to the initial simulated BCG signals. Here we consider two low-pass filters with a cutoff frequency of 25 and 40 Hz, respectively. Figure 6 already shows the impact that these filters have on and .
It has already been shown that the different time intervals measured on were only slightly affected by low-pass filters, down to a cutoff frequency of 25 Hz (Gómez-Clapers et al., 2013). Supplementary Figure 3 presents the Bland-Altman plots used to compare the results obtained with the scaling law, based on a low-pass filtered signal with a cutoff frequency at 40 versus 25 Hz. The mean differences observed on all the parameters studied are very close to 0, even though ̃ appears to be a bit overestimated with a cutoff frequency at 25 Hz, in comparison to 40 Hz. Even for this parameter, however, the mean difference is close to 10% of the mean value observed for ̃ in all the simulated cases. Such a difference can be considered relatively low in comparison with the wide range of values taken for the inverse of pulse wave velocity in the aorta (twice as large at 20 years old as at 80 years old). This confirms the potential of this scaling law to provide parameters that are only slightly affected by different body types and can allow the estimation of some vascular parameters such as pulse wave velocity in the aorta. Figure 3. Bland-Altman plots comparing results obtained based on a low-pass filtered signal with a cutoff frequency at 40 and 25 Hz for: a. ̃ ; b. ̃ ; c. ̃ ; and d. ̃. A positive difference indicates that the parameter of interest is larger with the cutoff frequency at 40 Hz than the one at 25 Hz. The mean difference is indicated by a solid line, while the confidence interval (mean +/-1.96 standard deviations) is indicated with dashed lines. The range on the vertical axis corresponds to +/-20% of the mean value observed for a given parameter. Colors indicate different simulated ages, while symbols indicate different stressed volumes. Several points can have the same color and symbol, when they correspond to different heart rates. The healthy young reference case is indicated by a black diamond.

7
Complementary results related to stroke volume

Determination of a scaling law linked to the integral of the kinetic energy
The scaling law giving 1 (Equation 34) has a physical basis described in the literature (Starr et al., 1939), and the one giving 2 (Equation 35) is an adaptation that may give more robust results, based on some in silico observations, as described in Section 3.4.2. The scaling law giving 4 (Equation 37) corresponds to another modification of Equation 34, based on the fact that sys should be less sensitive to inter-individual differences than the area under the I and J waves. However, it does not have a direct physical basis, as opposed to the scaling law giving 3 (Equation 36), which relies on the development of some relationships that are described hereafter.
Since the aorta is the largest contributor to , and thus kin , the following reasoning can be conducted based on Equation 22: with the blood flow rate through the aortic valve.
This means that: Furthermore, we can make the usual approximation that has the shape of a squared sinus of amplitude max and width . By integrating Equation A15 over the systolic phase (i.e., for kin , from the beginning of the cardiac cycle until the time corresponding to the end of the 'J wave), we find:

Effect of filtering on the results related to stroke volume
The metrics based on the amplitude of the BCG signals are also susceptible to be affected by filtering, and thus, as hypothesized, different body types. In this section, we are particularly interested about the estimates of stroke volumes obtained from and , according to the equations proposed in Section 3.4.2.
Supplementary Figure 4 presents the Bland-Altman plots used to compare the stroke volume estimates, based on an initial low-pass filtered signal with a cutoff frequency at 40 versus 25 Hz. The mean differences observed on all the stroke volume estimates are very small (below 10% of the mean value observed for all the simulated cases) and always positive. This indicates that the stroke volumes computed with a cutoff frequency at 25 Hz are underestimated, in comparison to 40 Hz. The scatter around the mean difference is very small for the four different formulas of stroke volume, but 3 appears to perform better, with a mean difference very close to 0. More generally, estimates based on seem to be less affected by filtering than those based on . This was to be expected, because the signal contains higher frequencies and the integration that leads to (and thus ) already plays a role of low-pass filter.
This short analysis confirms the potential of the suggested formulas to evaluate stroke volume and its changes with only a low impact of different body types on the estimates. Since the bias induced by low-pass filtering appears to be very stable, one could imagine that empirical laws could provide a correcting factor to apply to the BCG estimates, based on physical parameters of the subject.  37). A positive difference indicates that the parameter of interest is larger with the cutoff frequency at 40 Hz than the one at 25 Hz. The mean difference is indicated by a solid line, while the confidence interval (mean +/-1.96 standard deviations) is indicated with dashed lines. The range on the vertical axis corresponds to +/-20% of the mean value observed for a given parameter. Colors indicate different simulated ages, while symbols indicate different stressed volumes. Several points can have the same color and symbol, when they correspond to different heart rates. The healthy young reference case is indicated by a black diamond.

Effect of peripheral resistance on the estimation of stroke volume
Since peripheral resistance has an influence on stroke volume, we evaluate its impact on the quality of the estimates provided by the four suggested formulas. This analysis is conducted on the young healthy reference case, by gradually increasing the value of each systemic peripheral resistance (computed as 0 + 1 ) up to three times its initial value. The corresponding results are provided in Supplementary ). All the results are obtained after low-pass filtering of the acceleration signal at 25 Hz. R 2 : correlation coefficient computed on the whole dataset While 1 and 2 are visibly not able to correctly evaluate intra-individual changes of stroke volume caused by an increase in systemic peripheral resistance, 3 and 4 seem to keep a linear relationship with regard to stroke volume (R 2 = 0.95). However, the estimated changes are relatively underestimated in both cases (slope of 4.33 for 3 and 2.27 for 4 ).