{'fraction_non_alphanumeric': 0.02693922431027589, 'fraction_numerical': 0.015143942423030788, 'mean_word_length': 4.895403653506188, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic simulation of aortic valve stenosis using a lumped parameter cardiovascular system model with flow regime dependent valve pressure loss characteristics

Valvular heart diseases are growing concern in impoverished parts of the world, such as Southern-Africa, claiming more than 31 % of total deaths related to cardiovascular diseases. The ability to model the effects of regurgitant and obstructive lesions on the valve body can assist clinicians in preparing personalised treatments. In the present work, a multi-compartment lumped parameter model of the human cardiovascular system is developed, with a newly proposed valve modelling approach which accounts for geometry and flow regime dependent pressure drops along with the valve cusp motion. The model is applied to study various degrees of aortic stenosis using typical human cardiovascular parameters. The results generated with the proposed model, are compared to predictions using previously published valve modelling approaches and both sets of results are compared to typical local and global physiological parameters found in literature such left-ventricular systolic pressures, peak and mean aortic valve pressure drops and vena contracta velocities. The results show that the previously published valve models under predicts expected severely stenosed peak and mean transvalvular pressure drops by approximately 47% and 30% respectively, whereas the newly proposed model under predicts the peak pressure drop by 20% and over predicts mean pressure drop by 7%.


Introduction
Cardiovascular diseases (CVDs) are a growing social and healthcare burden on the world. Recent studies show that globally an estimated 19 million people died [1] from CVDs in 2020 compared to 17.9 million in 2017 [2]. The leading cause of CVD related deaths in high-income countries such as in Europe and Northern America is ischaemic heart disease (IHD), which accounts for approximately 50% of deaths from CVDs. In impoverished countries, such as those located in Africa, IHD accounts for less than 10% of deaths from CVDs. In Sub-Saharan Africa, valvular heart diseases such as rheumatic heart disease (RHD), bicuspid aortic valve disease and degenerative aortic valve disease, accounts for approximately 31% of CVD related deaths [3], [4]. Recent statistical projections also show that valvular diseases are set to double in the United States and Europe before 2050 [1].
Valvular heart disease is broadly categorised as obstructive and regurgitant lesions of the valve body. Valve malfunction, whether due to regurgitation or obstruction is associated with abnormal intracardiac haemodynamic behaviour.
Left untreated this can lead to reduced cardiac output or high ventricular pressures [5]. The ability to quantify the influence of valvular diseases on the local and global haemodynamics of the cardiovascular system would assist clinicians in evaluating patient specific cardiovascular performance and preparing personalised treatments [6]. One approach to non-invasively predict and study valve disease effects on local and global haemodynamic functions is to use electricalhydraulic analogue lumped-parameter mathematical models (LPMs) of the cardiovascular system [7]. LPMs can be used to quickly simulate variations in blood pressures, flow rates and volumes of various compartments that make up the circulatory system as a function of time.
[10], [11], [12], [13]. Garber et al. [6] stated that the heart valve modelling approaches predominantly used in literature, are merely simplifications of the actual valve mechanics and does not take into account the local haemodynamics such as mechanical forces and valve motion. Some researchers have developed more advanced valve models for LPMs, which included valve motion and areadependent pressure losses. Korakianitis and Shi [14], [15] studied the effect of aortic regurgitation and mitral stenosis using a LPM and an advanced valve model which accounted for valve motion and variable area effects on pressure losses. Similarly, Mynard et al. [16] proposed an advanced valve model which incorporates the valve motion as a function of the pressure forces. The proposed model used tuned coefficients to accurately capture the valve motion and pressure losses through the valves.
In the present work, a 0D valve model is proposed which accounts for valve motion and Reynolds number dependent pressure losses. The model, similarly to [15], solves the valve cusp dynamics as a set of differential algebraic equations (DAEs) and uses the valve position to calculate the pressure losses. To calculate the pressure drop through the valve, the domain is approximated as a gradual contraction (nozzle) rather than a typical orifice and includes the entrance, exit, frictional and local loss effects. The proposed model is applied to investigate the influence of increasing degrees of aortic stenosis on a cardiovascular system using a multi-compartment LPM. To investigate if the proposed level of valve modelling detail is necessary, the model outputs are compared to results generated using the approaches presented by [15]. According to the best knowledge of the authors the present work is the first instance of using a dynamic valve model with Reynolds number dependent pressure losses to study the local and global effects of aortic stenosis on a LPM of the circulatory system.
The computer models of the cardiovascular system, were developed using the Julia v.1.7.0 [17] programming language and the DifferentialEquations.jl [18] and NLSolve.jl open-source libraries.

Material and methods
In the current section, the mathematical models for the different components in the cardiovascular system seen in figure 1 will be discussed. The LPM developed in the present work consists of the four heart chambers and their respective downstream valves, the systemic circulation network and the pulmonary circulation network.

Modelling of heart chambers
The heart chambers act as non-linear pressure sources in the hydraulic network, which pressurises the low-pressure blood coming from the upstream veins or chambers. To simulate the pressure-volume-time relationships of the left ventricle (LV), left atrium (LA), right ventricle (RV) and right atrium (RA) the time-varying elastance model of Suga et al. [19] is used. The instantaneous static pressure inside the i th chamber (ventricles or atriums) is calculated using the relation shown in equation 1. The instantaneous volume V i (t) [mL] of the i th chamber is calculated using a simple mass balance over the respective heart chamber, as shown in equation 2.
In equations 1, 2: P i,0 [mmHg] is the unstressed chamber internal pressure, is the unstressed chamber volume, Q i,in (t) [mL/s] is the inlet flow rate into the chamber and Q i,ex (t) the exit flow rate out of the chamber. e i (t) is the time-varying elastance of the i th chamber and is estimated differently for the atriums and ventricles. For the ventricles the time-varying elastances are modelled using equation 3 [15] and the ventricular activation functions with equation 4 [10].
In equation 4, T [s] is the heart beat period, T 1 is the time at end of systole and T 2 is the time at end of ventricular relaxation. The time-varying elastances for the atriums are calculated using equation 5 and the atrial activation functions are estimated using equation 6.
In equations 5, 6: E i,max , E i,min are the maximum and minimum pressurevolume relations inside the atriums, T a is the time at onset of atrial contraction and D is the time of atrial relaxation. As seen above, in the time-varying elastance model equations, various empirical values are required, such as T 1 , Table 1 contains the empirical values used for the ventricles and atriums. In the current work the heart beat period is fixed at T = 1 [s].

Systemic and pulmonary networks
The systemic and pulmonary vasculature networks are modelled using five components each, as seen in figure 1. These components are the sinuses, arteries, arterioles, capillaries and veins [14]. Depending on the type of vasculature, the local haemodynamics is simulated using a combination of hydraulic resistance (R), blood inertia (inductance, L) and vessel compliance (capacitance, C) to capture the time-dependent pressure-volume behaviour of the specific vessel.
The flow rate through the aortic sinus is simulated using the ordinary differential equation (ODE) shown in equation 7. The time-dependent blood pressure inside the aortic sinus can be evaluated by solving the flow rate differential equation and equation 8 simultaneously [20], [14].
In the equations above, L AS mmHg · s 2 /mL is the blood inertia at the For the pulmonary sinus the flow rate and blood pressure can be simulated by solving equations 9, 10. The descriptions and units of the variables are similar to that of the aortic sinus above.
The systemic arteries are modelled similarly to the aortic sinus. The governing differential equations for flow rate and inlet blood pressure for the systemic arteries is given in equations 11, 12.
The variations in pulmonary arteries outlet flow rate and inlet blood pressure are simulated by solving equations 13, 14.
The systemic arteriole and capillaries are modelled using only hydraulic resistance components due to the rigidity and small diameter of the vessel walls [13]. The flow though these vasculature sections are assumed to be steady, The pressure drop through these sections are calculated using equations 15, 16.
Similarly, for the pulmonary network arteriole and capillaries the pressure drops are calculated using equations 17, 18.
The systemic and pulmonary veins are simulated as compliance-resistance components with negligible blood inertia. The governing differential equations for the systemic (equation 19) and pulmonary (equation 20) vasculature inlet pressures are shown below.
To determine the flow rate through the venous system the hydraulic resistance expressions are used. The driving pressure for the systemic vein flow rate is the pressure difference between the vein inlet and the right atrium chamber pressure [20]. For the pulmonary network, the driving pressure is between the inlet to the veins and the left atrium chamber pressure.
All parameters used in the vasculature model are shown in table 2. In the next section the modelling of the different valve modelling approaches will be discussed.

Heart valve models
The purpose of the four heart valves is to prevent reverse flow or regurgitation of blood back into the heart chambers, similar to the operation of check valves used in the process industry. The opening and closing processes of the valve cusps are governed by a combination of pressure, vortex and frictional forces [14] acting on the surfaces of the valves. In typical lumped parameter modelling of cardiovascular systems, heart valves are approximated as resistance components with a single diode component [22]. Ref.

Valve model 1
The simplified heart valve model presented by [14], estimates the flow rate of blood through the i th valve using equation 23, which is based on a typical orifice pressure drop relation [23]. In equation 23, A base is the base area of the inlet conduit to the valve, ρ l is the blood density, K is the flow coefficient, A r is the valve area opening ratio defined as  [14].
The valve area opening ratio A r shown in equation 23, is estimated using simple diode behaviour, where A r either takes a value of 0 or 1 depending on the pressure gradients across the specific heart valve, as shown in equation 25.
In equation 25, P in (t) [mmHg] is the instantaneous valve upstream static blood pressure and P ex (t) [mmHg] is the instantaneous valve downstream static blood pressure.

Valve model 2
The advanced heart valve modelling methodology presented by [14], also To find the time-dependent θ v (t) [rad/degrees] for each heart valve, four differential equations governing the valve cusp motion are solved. These dynamic equations, takes into account the pressure forces acting on the valve cusps, frictional forces due to neighbouring tissue resistance and the fluid-structure interactions of the vortexes forming downstream of the valves. Therefore, to simulate the rotational cusp motion of the i th heart valve, equation 27 is numerically solved.
In the equation above, F p,i (t) is the pressure forces, F f,i (t) the tissue frictional forces and F v,i (t) the vortex forces acting on the i th heart valve at time step t. I o,i is the mass moment of inertia of the heart valve cusps.
As discussed in the works of Korakianitis and Shi [14], [15] and Mynard et al. [16], the pressure force is the dominant factor which dictate valvular motion. Therefore, in the current work only the pressure forces are considered and frictional and vortex forces ignored (F f = F v = 0). To estimate the pressure forces on the cusps of the i th heart valve the constitutive relationship shown in equation 28 is utilised.
In the equation above, the valvular force coefficient K p rad s 2 ·mmHg is assumed to have a constant value of 5500 for all four heart valves. This force coefficient value was determined through manual model tuning to produce valvular motion similar to that found in literature. From the above discussion it is apparent that the proposed model, does not account for all the mechanical forces acting on the AV valves such as cord tensions. Nonetheless it has been shown that the valve model 2 dynamics adequately captures the motion of the AV valves.

Valve model 3
The challenge associated with implementing the valve modelling methodologies discussed in sections 2.3.1 and 2.3.2 in patient-specific LPMs, is that the flow (CQ) and valvular force (K p ) coefficients are manually tuned and, therefore, not applicable to a wide range of valve geometric and haemodynamic parameters.
In the present section, valve pressure loss and motion models will be presented that are based on basic valve parameters such as valve cusp thickness, cusp height, valve opening angle and instantaneous valve flow rate. Furthermore, the pressure loss calculations will distinguish between AV and semilunar valves due to the different configurations relative to the adjacent heart chambers. The   dθ The mass moment of inertias for the cusps are calculated using equation If one studies the geometry of the heart valve shown in figure 2 In the equation above K t is the total flow coefficient for the semilunar valves and consists of K c,SL the local contraction loss coefficient, K f,SL the semilunar valve friction loss coefficient and K se,SL the semilunar valve sudden expansion coefficient. The contraction loss coefficient is calculated as shown in equation 33 [27], [28].
In equation 33, α i (t) = π 2 − θ v,i (t) and β i (t) = A r (t) for the i th valve. The contraction ratio is calculated as λ i (t) = 1+0.622· αi(t) The frictional loss coefficient is calculated using equation 34.
The friction factor f i per semilunar valve is calculated using the well-known smooth surface Colebrook-White equation [23] shown equation 35. The Reynolds number used in the calculation of the friction factor is determined using equation 36, and is based on the velocity at the vena contracta of the semilunar valve.
The sudden expansion loss coefficient for the semilunar valves is calculated using equation 37 [28], where

Case studies
In the current work, the valve models presented in sections 2. For valve model 1, the area opening ratio A r is simply limited to a specified value and for valve models 2 and 3 the maximum opening angles are limited θ v,max,AO . Table 4, contains the limits applied to the aortic valve for the 5 case studies.

Numerics
To simulate the dynamics of the cardiovascular system, the differential- are solved simultaneously using an ODE system of equations integrator. For the current work the Tsitouras 5/4 Runge-Kutta method [34] is used with the relative and absolute tolerances set to 1 · 10 −4 and 1 · 10 −6 respectively. The maximum integrator iteration count of 10 6 was also specified.  To numerically integrate and solve the mentioned DAEs, certain constraints and conditions had to be incorporated such as valve motion limits, implicit equation solving and initial conditions. To incorporate the discontinuities arising from the valve motion limits the following condition is included in the simulation procedure of each valve

Results and discussions
To demonstrate that the proposed LPM, along with the three heart valve modelling extensions, can reproduce typical physiological haemodynamic parameters the results for case study 1 was compared to data from literature as seen in table 5. Note that M1, M2 and M3 corresponds to the LPM implemented with valve models 1, 2 and 3 respectively. The results in table 5 shows that the three models can capture typical physiological parameters of the human cardiovascular system. It is seen that M2 and M3 produce very similar results for case study 1, whereas M1 generates results with a higher cardiac output and higher blood pressures across the network.
These differences between models are due to the different valvular pressure drops calculated by each approach and are reflected in table 6.  and Shi [14]. The percentage differences between M1 and M2 LVSPs over the range of stenosis cases analysed stay relatively constant at 7.6%. M1 predicts on-average 9% higher stroke volumes compared to the values predicted by M2.
Therefore, M2 predicts lower cardiac output, which highlights the importance of including valve motion effects seeing as cardiac output is a critical clinical parameter [8]. The LVSP and stroke volume trends observed in figure 5, indicates that the M3 modelling approach predicts higher pressure drops (ventricular afterload) in the stenosed valve compared to the M1 and M2 models.    proposed by [14], which results in M3 calculating higher pressure drops through the mitral valve when compared to the values of M2 (table 6). Studying the aortic valve data, it is observed for case study 1, that the calculated flow coefficients by M3 are significantly larger than the constant value used by M1 and M2. This explains why in table 6, the M3 model predicts a lower aortic valvular pressure drop compared to M2. In figure 8 (b), for the moderate stenosis case we seen that the flow coefficients calculated by the M3 approach, are lower compared to the constant value used in M1 and M2. The lower flow coefficient values result in larger pressure drops over the valve during valvular ejection.
For the very severe aortic stenosis case it is observed that the flow coefficients are lower than in the moderate stenosis case. This reduction in flow coefficients as the degree of stenosis is increased, is due to the higher sustained blood flow Reynolds numbers in the valve which leads to larger frictional and geometrical contraction losses.
The approximate maximum and minimum flow coefficients calculated by the M3 modelling approach are 650 and 150 respectively, over the range of case studies. This gives an average flow coefficient value of 400 which is close to the tuned value proposed by Korakianitis and Shi [14] (350 and 400). The advantage of using the M2 approach, is that it solves significantly faster than the M3 approach. The M1 and M2 modelling approaches, simulates 10 s, in approximately 15 s, whereas M3 takes 10-15 minutes.

Conclusions
The results in the present work shows that the newly proposed valve modelling approach, predicts higher pressure drops at more severe degrees of aortic stenosis when compared to the results of the models from literature. This is due to the fact that the models from literature are insensitive to changes in blood flow Reynolds number through the valve. Comparing the calculated ranges of pressure drops over stenosed valves with experimental data from literature, shows that the proposed model captures the expected pressure drops more accurately compared to the models from literature.