The time-varying elastance of each cardiac chamber, e ( t ) (Equation 1), is described using a double-Hill function (Stergiopulos et al., 1996; Mynard et al., 2012; Broomé et al., 2013) that captures both the contractile and passive mechanical properties of the myocardium: e t = k g 1 1 + g 1 1 1 + g 2 + e min v , (1) where g 1 = t α 1 T m 1 , g 2 = t α 2 T m 2 . (2) In Equation 2, T denotes the duration of the cardiac cycle. The dimensionless parameters α 1 , m 1 , α 2 , and m 2 define the shape of the elastance curve during the contraction and relaxation phases, respectively. The scaling factor k is chosen such that the maximum elastance satisfies max ( e ( t ) ) = E max : k = e max v e d , q − e min v max g 1 1 + g 1 1 1 + g 2 . (3)

In Equation 3, the maximum elastance, e max , is modulated by the end-diastolic volume v e d and chamber-specific output flow q , thereby incorporating the Frank–Starling mechanism (Sun et al., 1995) as shown in Equation 4: e max v e d , q = E max ⋅ 1 − v e d V e d , max 4 ⋅ 1 − q Q max . (4) The minimum elastance, e min ( v ) (Equation 5), describes the passive exponential pressure-volume relationship during diastolic filling (Chung et al., 1997): e min v = E min ⋅ e σ ⋅ v − v 0 . (5) Parameters related to the cardiac function are reported for each chamber in Table 1.

Interventricular interaction is modeled through a septal elastance, E s v , which couples the pressures and volumes of the left and right ventricles (Sun et al., 1997). The left ventricular pressure, p l v (Equation 6), is expressed as the sum of its effective elastance contribution and the cross-talk pressure originating from the right ventricle: p l v = E s v E s v + e l v ⋅ e l v ⋅ v l v − V 0 , l v + e l v E s v + e l v ⋅ p r v , (6) where E s v denotes the septal wall elastance, e l v is the left ventricular elastance as defined in Equation 1, and p r v is the right ventricular pressure. This relation captures septal shifting during the cardiac cycle and its impact on the effective pressure within the left ventricle. The right ventricular pressure contribution and atrial septal interactions are modeled in a similar way.

The pericardium is modeled using an exponential pressure–volume relationship (Equation 7) to represent the mechanical constraint it imposes on cardiac filling: p pc v = p p c , 0 + K p c ⋅ e v p c − V p c , 0 φ p c , (7) where p p c , 0 denotes the minimum pericardial pressure, K p c is the pericardial pressure constant, V p c , 0 represents a volume offset, φ p c is the pericardial volume constant, and v p c corresponds to the total pericardial volume, defined as the sum of the volumes of the cardiac chambers and the pericardial fluid volume V p e .

The pressure drop across a cardiac valve is governed by Equation 8, a nonlinear inertial-resistive model (Mynard et al., 2012): Δ P = B ⋅ q ⋅ | q | + L ⋅ d q d t , (8) where the resistive and inertial coefficients are defined by Equation 9: B = ρ 2 A eff 2 , L = ρ l eff A eff , (9) where ρ = 1,060 g/ c m 3 is the blood density, and l eff is the inflow length, which is set to be equal to the instantaneous opening diameter of the valve (Broomé et al., 2013). The effective valve area A eff (Equation 10) changes dynamically based on a gating variable ζ ( t ) , which models the valve’s opening and closing dynamics: A eff t = A eff , max t − A eff , min t ζ t + A eff , min t . (10) The dynamics of ζ ( t ) are governed by Equation 11: d ζ d t = 1 − ζ K vo Δ p , if  Δ p > 0 ζ K vc Δ p , if  Δ p < 0 , (11) where K v o and K v c are the rate coefficients for valve opening and closing, respectively. Cardiac Valves parameters are reported in Table 2.

Vascular dynamics of the venous and arterial compartments are modeled using a four-element Windkessel model. The flow update at the generic i t h vascular compartment is computed as in Equation 12: d q i d t = p i + 1   − q i  ⋅ R − p C − p i t − q i - 1 − q i ⋅ Ω + Δ p h i / L , (12) where q i - 1 and p i + 1 i + 1 denote the intravascular blood flow and transmural pressure of the previous ( i - 1 ) and subsequent ( i + 1 ) compartment, respectively. R , Ω , and L are the resistance, viscoelastance, and inertance of the blood vessel, while Δ p h i is the gravitational hydrostatic contribution. The transmural pressure is a sum of the viscoelastic contribution q i - 1 − q i ⋅ Ω , the intra-thoracic pressure p i t and the pressure p C (Equation 13) which follows an exponential volume-pressure relation (Dahn et al., 1970): p C = P 0 e v − V 0 φ , (13) where V 0 is the volume at mean pressure P 0 , while ϕ is a vessel-specific parameter which determines the non-linearity of the pressure-volume relation. The viscoelastic damping term Ω (Equation 14) depends on the vessel’s and fluid inertia: Ω = λ E L , (14) where λ = 0.5 is an additional damping factor that helps in maintaining numerical stability throughout the simulation.

Geometric and vessel material properties define the resistance (Equation 15), inertia (Equation 16), and elastance (Equation 17) of each compartment: R 0 = 8 ⋅ η ⋅ l π ⋅ r 0 4 ⋅ n , (15) I 0 = ρ ⋅ l π ⋅ r 0 2 ⋅ n , (16) E 0 = Y inc  ⋅ h 2 ⋅ π ⋅ r 0 3 ⋅ l ⋅ n , (17) where l , r and h are the length, radius and wall thickness, η = 0.00024 mmHg ⋅ s is the blood viscosity, Y inc the Young’s modulus, and n is the number of parallel vessels, lumped into each compartment. The actual segmental elastances, resistances and inertances are updated in each calculation step based on the volume and radius, assuming constant vessel length and thickness.Vascular compartment - specific paraemters are reported in Table 3.

The time-averaged aortic pressure p ̄ ao is regulated via a baroreflex mechanism (Fois et al., 2022), modeled using sigmoidal activation functions for the sympathetic ( n s ) and parasympathetic ( n p ) neural firing rates as shown in Equation 18: n s p ̄ ao = 1 1 + p ̄ ao p ̄ ao , t g ν , n p p ̄ ao = 1 1 + p ̄ ao p ̄ ao , t g − ν , (18) where ν = 7 defines the steepness of the neural response, and p ̄ ao , t g is the target aortic pressure, set to 100 mmHg to match the mean arterial pressure (MAP) of the model in the supine position. These neural signals modulate a generic cardiovascular variable y m ( t ) according to Equation 19: d y m d t = 1 τ m − y m + α m n s p ̄ ao − β m n p p ̄ ao + γ m , (19) where α m , β m , γ m , and τ m are parameters characterizing the baroreflex dynamics specific to the regulated variable. Controlled variables include HR, systemic vascular resistance (SVR), left ventricular maximum elastance E l v , m a x , systemic venous compliance C sys,ven , and systemic unstressed venous volume V u n , s y s , v e n .

Cardiopulmonary regulation of the time-averaged right atrial pressure p ̄ ra is implemented analogously, with a target pressure p ̄ ra , t g = 5 mmHg. This mechanism modulates pulmonary vascular resistance (PVR), right ventricular maximum elastance E r v , m a x , pulmonary venous compliance C pul,ven , and pulmonary unstressed venous volume V u n , p u l , v e n .

Baroreflex gains were empirically selected to match the cardiovascular blood pressure response as reported in previous work (Fois et al., 2022), which is based on in-vivo data (Coonan and Hope, 1983; Smith et al., 1994; Blomqvist and Stone, 1991). The relative contributions of HR, SV, CO, TPR and LVemax were also controlled via the baroreflex gain. Parameters related to the baroreflex model are reported in Table 4.

In walking, two accelerations act on the vascular system: gravity and body acceleration. The effect of hydrostatic pressure induced by gravity and the hydrodynamic pressure induced by body acceleration is modeled as an additional pressure head in each vascular segment: Δ p h i = ρ ⋅ l i ⋅ a ⋅ cos α t , (20) where ρ is the blood density, l i is the length of the i -th segment, a is the vertical acceleration acting on the i -th segment, and α ( t ) is the time-dependent angle between the body axis and the gravity vector. The length l i of each vascular compartment is assumed to match the corresponding values reported in Table 3.

Intra-thoracic pressure p i t (Equation 22) is also adjusted in the standing position, transitioning linearly from − 2.5  mmHg in the supine position to − 6.5  mmHg in the upright position, as a function of posture angle: p i t α = p sup i t + p up i t − p sup i t ⋅ α 90 ° (22)

where p sup i t = − 2.5  mmHg and p up i t = − 6.5  mmHg .

We conducted a tilt test simulation to validate our integrated baroreflex and gravitational model against published physiological data (Fois et al., 2022). For this purpose, mean values of main pressures, including MAP, systolic arterial pressure (SAP), diastolic arteral pressure (DAP) and central venous pressure (CVP) were compared against mean values from experimental literature and computational data from Fois et al. (2022). Percentage change in HR, SV, CO, total peripheral resistance (TPR), and E l v , m a x was also compared.

Postural transition was modeled as a 90° head-up tilt from the supine position, following the methodology described in Heldt (2004).

The postural angle α ( t ) (Equation 23) was modeled using a smooth cosine transition to emulate the physiological tilt maneuver: α t = α sup , t < t tilt , 0 α sup + α up − α sup 2 1 − cos π t − t tilt , 0 Δ t tilt , t tilt , 0 ≤ t ≤ t tilt , 0 + Δ t tilt α up , t > t tilt , 0 + Δ t tilt (23)

where α sup = 0 ° is the supine angle, α up = 90 ° is the standing angle, t tilt , 0 is the tilt initiation time, and Δ t tilt is the duration of the postural change.

Synthetic acceleration signals were computed by repeating one vertical acceleration cycle and were provided as input to the model. In all simulations, the HR was fixed at 70 bpm. The pressure in the descending aorta was selected for this analysis.

Initially, the SR was set equal to the HR, with five distinct phase shifts to the cardiac cycle. The first phase shift occurred when the peak acceleration coincided with mitral valve closure, the third phase shift corresponded to the alignment of peak acceleration with ventricular filling, marked by maximal mitral valve opening, and the last phase shift when the peak acceleration occurred just before the mitral valve closure of the subsequent heart beat. The second and fourth shifts occurred during the transitions between systole and diastole and diastole and systole.

Subsequently, the SR was increased beyond HR by 10%, 20%, and 30%. The difference in frequency, naturally iterates through the different phase shifts of the two physiological systems. Fast Fourier Transform (FFT) was performed to compare the frequency content of the pressure signals between the standing and walking conditions. FFT resolution was set as 0.02 Hz.

The experimentally measured chest acceleration and HR data from the subjects were used as inputs to the computational model. The simulated pressure in the descending aorta was then compared to the brachial arterial pressure.

To determine the phase synchronization between step events and the cardiac cycle, the phase of each step was calculated as the percentage of the RR interval at which the heel strike occurred. Heel strike timing was identified from insole data as the point at which average pressure reached a certain threshold, calculated with Cometa proprietary software. Cardiac contraction was defined as the R-peak in the ECG for experimental data and as the mitral valve closure in the simulated data. A phase of 0% corresponds to a step occurring simultaneously with cardiac contraction, while 100% indicates a step coinciding with the subsequent contraction.

Pulse waveforms were segmented to each phase, starting from heart contraction to 35%RR after the subsequent contraction. The mean waveform for each phase was then computed for further analysis. The experimental waveform was shifted 670 samples for alignment purposes.

At rate-equivalence, inertial effects associated with walking are evident in the normalized pressure waveforms of the descending aorta, as illustrated in Figure 3. The Figure shows increasing phase shifts (arrows) between the heart contraction (red) and the maximum body acceleration (blue). Compared to the baseline waveform while standing (Figure 3A), walking introduces an additional hydrodynamic pressure component (Figure 3B), resembling the waveform shape of the chest acceleration. As the acceleration peak shifts with respect to the cardiac contraction, the resulting inertia-induced pressure augmentation also shifts towards the diastolic part of the waveform to the right. Depending on the phase shift between cardiovascular and locomotor activity, and thus on the timing of this interaction, the increase is observed either in systolic or diastolic pressure, (Figure 3C).

When the SR exceeds the HR, the phase is no longer constant and we expect higher order system interactions between two systems of similar frequency. In the time-series signal, we observe a low frequency component in the pressure waveform, characterized by a frequency lower than both HR and SR, Figure 4A. The MAP also reflects this interaction, exhibiting a sinusoidal pattern, with low frequency.

The frequency spectrum in Figure 4B reveals distinct peaks at the SR the HR and their higher order harmonics. As presented in Table 5, we can link most frequency peaks to the imposed HR and SR. The dominant low frequency is associated with the lower frequency observed in the MAP. The value is equal to the difference between SR and HR across all conditions tested, which suggests that the signal results from a ‘beating’ interaction between the HR and SR. However, the addition of the two frequencies does not result in a frequency peak. Additional low-frequency components emerge in the range between the ‘beating’ frequency and the HR. The low magnitude of these components (f < 0.1 mmHg) suggests that they might be associated with harmonics of the beating frequency or could result from numerical aliasing effects.

A total of 1.120 waveforms for P1 and 750 for P2 were included in this analysis. The similarity index and cumulative distribution function are shown in Figure 8 for both participants. The morphology results indicate an increase in similarity index when hydrodynamic pressure is incorporated in the model for both subjects (K-stat 0.123 vs. 0.029, P1 and 0.164 vs. 0.059, P2). The distributions are closely aligned for the two extremes 0 and 1, representing pressure onset and systolic peaks. Larger difference is instead observed for pressure values between 0.4 and 0.6, most likely to occur at the diastolic phase of cardiac cycle. It is further noted that all pairwise comparisons between simulation and experiment still indicate that the distributions remain significantly different (p < 0.05).