First blood: An efficient, hybrid one‐ and zero‐dimensional, modular hemodynamic solver

Low‐dimensional (1D or 0D) models can describe the whole human blood circulation, for example, 1D distributed parameter model for the arterial network and 0D concentrated models for the heart or other organs. This paper presents a combined 1D‐0D solver, called first_blood, that solves the governing equations of fluid dynamics to model low‐dimensional hemodynamic effects. An extended method of characteristics is applied here to solve the momentum, and mass conservation equations and the viscoelastic wall model equation, mimicking the material properties of arterial walls. The heart and the peripheral lumped models are solved with a general zero‐dimensional (0D) nonlinear solver. The model topology can be modular, that is, first_blood can solve any 1D‐0D hemodynamic model. To demonstrate the applicability of first_blood, the human arterial system, the heart and the peripherals are modelled using the solver. The simulation time of a heartbeat takes around 2 s, that is, first_blood requires only twice the real‐time for the simulation using an average PC, which highlights the computational efficiency. The source code is available on GitHub, that is, it is open source. The model parameters are based on the literature suggestions and on the validation of output data to obtain physiologically relevant results.

Recent years brought numerous solvers to model arterial flow in one dimension (1D), some in zero dimension (0D) too; see Table 1. The number of artery vessels is typically between 50 and 100. Although most of the considered papers also investigate the cerebral arteries in the brain, only nine out of 25 models the heart and seven the coronary arteries. While a lumped (or 0D) model can mimic the essence of the human heart with basic linear elements, for example resistor or capacitor, the coronary arteries are typically modelled in 1D and 0D combined. 44 However, the literature still lacks analysing the behaviour of the venous system. The model in this paper consists of the systemic and cerebral arteries with 99 vessels overall, and lumped models mimicking the heart and the peripherals.
In the last decade, two major articles presented comprehensive studies about different numerical solvers for onedimensional hemodynamics. Wang et al 45 investigated four schemes (MacCormack, upwind scheme, Taylor-Galerkin, and discontinuous Galerkin), and all of them are capable of predicting pulse wave propagation in physiologically normally conditioned humans (i.e., no shock waves). However, the computational times differ significantly, favouring the Taylor-Galerkin and the upwind scheme. Boileau et al 46  trapezium rule method. Their conclusions coincide, that is, every technique can resolve the pulse waves with appropriate accuracy. As Table 1 shows, the finite element (FE) and the finite difference (FD) methods dominate in 1D. Even though finite volume (FV) techniques are widespread in higher dimensions, they are barely visible in low dimensions. The most commonly used discretisation approach is a version of the Galerkin method, the original Taylor-Galerkin, or the discontinuous Galerkin. The simple first-order implicit method or the second-order Lax-Wendroff (or MacCormack) technique is applied for FD. Although the method of characteristics (MoC) can have some advantages over the former ones (e.g., computational time, handling shock waves), similarly to the FV, its visibility in the literature is negligible. As later described in detail, the MoC is used in this paper.
This paper aims at presenting computationally efficient and accurate numerical methods of our 1D-0D combined solver, capable of modelling human hemodynamics. The description focuses on the numerical issues and challenges in connection with the solver, with particular emphasis on how the combination of a 1D MoC solver and a nonlinear 0D model works properly with high stability and accuracy. We would also like to highlight the computational time of the MoC solver compared to popular FE and FD solvers. The key features of the presented solver, first_blood, are • open source, 47 • nonlinear 0D elements, and • modular in terms of topology both in 1D and 0D.
The model building for the human arterial system is also presented, emphasising the difficulties with parameter adjustment.
The paper is organised as follows. Section 2 introduces the basics of the mathematical modelling of the onedimensional human hemodynamics, their numerical solution with the method of characteristics, the general lumped solver, and the relationship between the 1D and 0D solvers in detail. The model building and the parameter tuning are presented in Section 3. As this paper aims at introducing the mathematical background and the numerical challenges behind the simulations, the results obtained with the model and the discussion are only briefly presented in Section 4. Section 5 summarises the work.

| Method of characteristics
The traditional equations of 1D fluid dynamics are solved, that is, the conservation of mass and momentum. Besides, the arterial wall is viscoelastic, described by the Poynting-Thomson model. The mass balance equation with a varying cross-section area describes the connection between the velocity and the cross-section along a pipeline. 48 A represents the surface area of a cross-section, v the cross-section average velocity, t the time, and x the spatial coordinate. Assuming laminar flow and circular cross-section, the momentum equation 48 connects the velocity v and the pressure p. In Equation (2), ρ is the density, ν the kinematic viscosity, and d the diameter. The Poynting-Thomson Equation (3) describe the wall material behaviour, a viscoelastic model with two springs in series and a damper parallel to one of the springs, see Figure 1. σ represents the tangential mechanical stress, ε the radial strain, E the Young modulus, and η the damping coefficient. B ardossy et al 3 presented the complete solution of Equations (1)(2)(3)(4)(5) with the method of characteristics and explicit Euler techniques. This paper briefly summarises the main steps and assumptions.
1. The Poynting-Thomson model is solved using the definition of the radial stain ε ¼ dÀ d n ð Þ=d n , the circular crosssection A ¼ d 2 π=4, the thin-walled pressure vessel equation σ ¼ pd=2s (where s is the wall thickness), material continuity d n s n ¼ ds (assuming that the Poisson coefficient is approximately 0.5), and assuming ε is small (i.e., ε 2 ≈ 0). Here, s is the wall thickness and the subscript n denotes nominal (unloaded). 2. Using the solution of the material model, the continuity equation can be reformulated by eliminating the crosssection A, then the strain ε, and the only remaining field variables are the velocity and the pressure in a partial differential equation (PDE) system with only two equations. 3. The method of characteristics transforms the PDE into two ordinary differential equations (ODE), and finally, the ODEs are discretised with explicit Euler.
The final form of equations are linear and algebraic, where J represents the source term F I G U R E 1 Mechanical model of the Poynting Thomson body. 3 E 1 and E 2 are the Young moduli, η is the damping factor, σ the tangential mechanical stress, and ε the radial strain containing the laminar pressure loss and the viscoelastic effects. The time step is denoted by Δt. The two main unknown variables (pressure and velocity) are at level P, and X stands for either R or L, that is, right and left, depending on the characteristic equation. For notations, see Figure 2. The wave speed describes the velocity of propagation information (e.g., pressure waves) in a motionless fluid. The parameter β mimics the effect that higher pressure increases the elasticity of arteries. 48 Equations (6) and (7) are linear and algebraic; thus, there is only one unique solution for the velocity and the pressure. Due to the characteristic lines, the solutions (point P) do not intersect with the equidistant grid. Linear interpolation can overcome such an issue: first, in time back to the new time level t iþ1 , then spatially to the equidistant grid. Once the pressure is calculated, the formulae derived from the Poynting-Thomson model can determine the strains. For the first time step, initial conditions must be defined, which can be p 0 atmospheric for the pressure and 0 mL/s for the volumetric flow rate. The initial diameter can be nominal, and the strain zero everywhere. In a general time step, every field variable is known at the ith time level, and the goal is to calculate the variables at the next, i + 1th time level. The time step Δt is predefined by the characteristic lines, depending on the exact length of the vessel and the number of inner calculation points.
where i is of the inner division points and i 0, N ½ . This formula also ensures the fulfilment of the Courant-Friedrich-Lewy condition, that is, each intersection of the characteristic lines must be above the new time level in a space-time coordinate F I G U R E 2 Visualising the method of characteristics along a single vessel. Red lines indicate the left (L), blue lines the right (R) characteristic lines, and P is their intersection. Δt is the time step, Δx the grid size, point T is the spatial coordinate at the old time level, point S at the new time level, and t i , t iþ1 the new time level system (otherwise, information is lost), see Figure 2. Although it ensures the Courant number is one, the time step is different for each vessel due to the different Δx grid size and pulse wave velocities. On average, the time step is around 0.75 ms.

| Boundary conditions
Boundary conditions (BC) must be defined at both ends of an artery, as there is only one characteristic line there. Three types of BC can be applied: • inner boundaries in MoC model, • prescribed boundary values in time, • or boundary conditions with a lumped model.

| Inner boundaries in MoC model
Considering the firstone, that is a junction between two or more arteries, the continuity equation is applied, where q i is the volumetric flow rate in the ith vessel, δ i indicates the direction of the vessel, and l is the number of vessels joining in the node. The second term models the neglected smaller arteries, muscles, and blood supplies for other tissues by a leakage with a linear resistance coefficient R, moreover, p P is the nodal pressure, and p 0 is the atmospheric pressure. Although the surroundings might have an overpressure, it is neglected and modelled as atmospheric. Furthermore, the characteristic equations (Equation 6 or 7) must be solved l times, where the unknown velocity v P ¼ q P =A P , and the cross-section is also a function of the pressure p P , where C 1 and C 2 are constants at a given time step originating from the viscoelastic property, see Equations (3)- (5), The iteration of the cross-sectional area is critical, as neglecting it (e.g., using the value from the previous time step) could result in numerical instability and inaccuracy. Overall, there are 2l þ 1 nonlinear, algebraic equations in each node: l of Equation (16), l of Equation (6) (or Equation 7), and the continuity (Equation 15), while the unknowns are the nodal pressure p P , l of volumetric flow rates q P and l of cross-sectional areas A P . The solution can be obtained computationally efficiently by expressing the volumetric flow rate from Equation (6) (or Equation 7) and Equation (16) as a function of the junction pressure, then substituting it into the continuity. It results in a single nonlinear, algebraic equation for the nodal pressure. Newton's iteration converges to the solution for such an equation within a few iteration steps.

| Prescribed boundary values in time
The prescribed boundary value can model the patient-specific measurement for pressure, volumetric flow rate or velocity. These are typically applied as upstream BCs. If the velocity is prescribed in time, the pressure can be directly computed with the characteristic equation (Equation 6 or 7). As far as the pressure, or the volumetric flow rate is considered, both the characteristic equation (Equation 6 or 7) and the cross-section equation (Equation 16) must be solved simultaneously. Similarly to the junction, Newton's technique can solve such nonlinear, algebraic equations efficiently.

| Boundary conditions with a lumped model
As different models are required to mimic the behaviour of organs, such as the heart or the peripherals, a general lumped or zero-dimensional (0D) solver is practical. Table 2 shows the types of elements a 0D model might include. The components are linear or nonlinear, and either an algebraic equation or an ODE must be solved for each. For the solution, first, the implicit Euler discretisation scheme transforms the differential equations into algebraic ones. The implicit approach is crucial to preserve the numerical stability of the solver. In addition to a large number of equations (up to 20), the characteristic equation (Equation 6 or 7) and the change in cross-section (Equation 16) must be solved for each joining 1D pipe. The latter is necessary because neglecting the cross-section change would lead to inaccuracies and numerical instabilities. However, Newton's method can cope with such a nonlinear, algebraic system.
The elastance is a speciality for modelling the effect of heart contraction. The elastance is the function of time and describes the pressure-volume relationship of the left chamber. It is physiologically general, as it applies to every healthy adult regardless its age, weight or living habits. 49 Figure 3 depicts the curve shape as the function of time. There are only two parameters: the peak H max and the minimum H min . 50 It models the heart contraction physically as a capacitor with varying capacitance. The fluid dynamic meaning could be a tank with a fluctuating base area.
Although the elastance function seems mathematically convenient-continuous, smooth, and differentiable multiple times-numerical instabilities occur during the solution of its equation. The appearance of the dH t ð Þ=dt over H t ð Þ is the root cause of this issue as the minimum value of H t ð Þ is low, while the derivative can be high, especially in the beginning of the cycle (see Figure 3). The matrix of the Jacobian during Newton's iteration has a high conditional number, indicating numerical instability. The problem can be resolved by eliminating the original elastance equation and introducing a new variable: y ¼ Δp=E. 50 Instead of solving the original elastance equation from Table 2, a new equation  from Table 3 is solved, avoiding numerical instabilities. Although it keeps the conditional number low, two additional linear, algebraic equations must be solved relating the new variable y to the pressure p at the starting and ending node of the elastance element.

| Algorithm design
A typical hemodynamic model contains numerous vessels with different lengths, from a couple to hundreds of millimetres. Since the range of the time steps can be high, it is convenient for each edge to have its own individual time T A B L E 2 Elements of the 0D model with their continuous equations and parameters. R is the resistance, C the capacitance, H t ð Þ the elastance, L the inductance, and V a constant pressure difference, Δp the pressure difference, q the flow rate Notation Type Continuous Eq.

Diode open
Δp ¼ Rq and q > 0 D Diode closed q ¼ 0 and q < 0 vector. Otherwise, the interpolation back in time could ruin the computational efficiency. Figure 4 depicts the layout of a hemodynamic simulation. The detailed description is the following. Loading the data from files and filling up every field variable. The nominal values (e.g., for the diameter), zeros (e.g., for the pressure) or results from previous simulations, if available, can be used as initial conditions. 2. Equation (14) determines every new time step for all arteries in the model. 3. The lowest new time level must be found with the corresponding jth index of the edge. Only the chosen jth edge will be updated to its new time level. Typically in each cycle, a different vessel is chosen. 4. Two independent conditions can stop the simulation based on engineering decisions. The first one is the prescribed simulation time, and if every artery reaches it, the simulation stops. The second one is the dynamic monitoring of the periodicity of the simulation results. The root mean square difference between the upcoming cycles quantifies the goodness of the periodicity. Although the model has many temporal histories, checking the aorta pressure signal provides sufficient accuracy and does not generate significant additional computation. The user can pick the stopping condition, and if the chosen one is fulfilled, the simulation stops. 5. If the simulation keeps running, the pressure and velocity fields are updated in the jth edge, which is at the lowest new time level. First, the inner points are calculated using Equations (6) and (7), see Section 2.1. Second, the exact equations for the boundary conditions depend on the surrounding: other arteries, prescribed field variable, or 0D model; see Section 2.2 for details. 6. Since every newly calculated P point is at a slightly different time level, linear interpolation can adjust them back to the lowest new time level with an equidistant grid; see Figure 2. Every other field variable can be updated at the new time level with the material Equations (3)-(5) from the pressure. 7. Last, the new time steps at the jth artery are determined, and then the code returns to Step 3. Table 4 contains all the general parameters of the models. The topology of the arterial system is based on the system presented. 11 Besides the main arterial vessels, the brain vessels, particularly the Circle of Willis and the neighbouring arteries, are also considered in detail; see Figure 5 top side. Moreover, the spine and the arteries of the organs around the spine are also part of the model. The heart and the peripherals are built using 0D models, 11,51 see Figure 6.

| HUMAN ARTERIAL SYSTEM
The standard, three-element Windkessel model mimics the peripheral effects of the arterial system. The structure of the heart is more complex: the mitral and aortic valves are modelled by an ideal diode and a resistance in series. While the diode is a model for the check valve, that is, it does not allow any backflow, the resistor models the pressure loss at the valve. The two inductors mimic the inertia of the fluid. A battery or "voltage" component prescribes the left atrium pressure. Keeping the left atrium pressure constant in low-dimensional modelling is an acceptable approximation. The elastance, that is, a capacitor with variable capacitance, is responsible for modelling the contraction effect of the heart. Interestingly, it is the same function in time for all patients, regardless of age, weight or race.
The most considerable challenge for building a proper model to describe the arterial system is not the topology but determining the countless parameters. Each artery has five geometrical (length, proximal and distal diameter, proximal and distal thickness), three material wall parameters (two elasticities, damping factor), and the number of inner division points is the ninth parameter. Moreover, the peripherals are modeled with a three-element Windkessel model; three independent parameters are at each artery ending. Finally, the heart is built using two resistors, one capacitor, one elastance and two inductors; overall, it means seven independent parameters. There are 99 arteries and 44 peripherals, meaning the number of all parameters is 99*9 + 44*3 + 7 = 1030. The question is how to choose these to have physiologically relevant outputs. With the help of imaging tools (e.g., computed tomography or magnetic resonance imaging), physicians can measure the geometrical parameters with high accuracy. The data found in Reymond et al 11 are used in the current research. Although Blanco et al 26 presented an approximation of the wall thickness via a power law function of the radius, here it is set to 10% of the diameter for all vessels as a rough estimate according to preliminary studies. 3 The value of the elasticities is critical, as it determines the wave propagation speed (Equation 19), thus the pulse wave velocity (PWV). Neglecting the radial strain Equation (10) simplifies to a n ¼ ffiffiffiffiffiffiffiffiffi Recently a study presented measured PWV values as a function of the nominal diameter. 52 A fitted approximation was used to achieve the physiologically proper PWV values (see Figure 7). The equation mathematically is a power-law function, which was also proposed in the literature, 52 that is where the α i parameters are determined by the least squares method. First, Equation (19) and the definition of the pulse wave velocity determine the E 1 Young modulus. Second, the spring parallel to the damper, E 2 , is set to 1:8 times of E 1 according to previous studies. 3 Besides the parameters, the number of inner grid points for every vessel must be determined. The primary goal is to accurately approximate the governing fluid dynamic equations with the least possible calculation points to maintain computational efficiency. The secondary goal is to obtain a uniform time step distribution, ensuring uniformity of accuracy. The time step is influenced primarily by the length of the artery, the number of division points and the PWV. Since the length is prescribed and the PWV is estimated according to Figure 7, a chosen number of division points determines the time step. The short arteries in the brain make it difficult to achieve uniform distribution because-based on preliminary simulations-at least five division points per edge are necessary to keep the accuracy high. As the vessel length decreases, so does the spatial grid, which also reduces the available time step, see Equation (14). For total uniformity, numerous calculations would be required with tens of thousands of points,  while the accuracy would be unnecessarily high. However, once an acceptable relative distribution is achieved, the overall number of inner points is decreased in the longer edges until the approximation error is negligible. Finally, the whole arterial system contains 1101 division points after some manual tuning and engineering decisions. Figure 8 demonstrates the effect of double and roughly half the grid points showing slight differences (below one mmHg) in the aortic pressure. Determining the parameters of the lumped models, such as the heart or the peripherals, is critical and challenging. None of the parameters can be directly measured or estimated. The only suitable method here is tuning the parameters to have a physiologically reasonable result, for example aortic pressure or cardiac output. The exact values for this study come from Reymond et al 11 for the peripherals, while the heart parameters are from Ferreira et al. 49 Finally, after setting every parameter, the initial simulations yielded correct temporal history forms, that is, the shape of the field variables was physiologically sensible; however, the amplitudes were not. Some additional tuning was necessary to increase the physiological correctness of the model. The calibration affected the peripheral Windkessel models and the heart model. Since the purpose of the paper is to demonstrate the 1D-0D solver first_blood, the process only included manual tuning based on engineering decisions and preliminary simulations.
The ethical statement is not relevant since there were no human measurements, only simulations with virtual models.

| Computational efficiency and accuracy
Although the computational time of low-dimension is hardly mentioned in the literature, it is a crucial aspect of solvers. A comprehensive study was published with four discretisation schemes: Taylor-Galerkin, monotonic upwind scheme for conservation law, MacCormack, and local discontinuous Galerkin. 45 Moreover, recently, Ghitti et al 43 published an efficient solver. Table 5 summarises the main aspects of computational efficiency. The table does not represent a comprehensive study, as not only the numerical schemes are different, but the environment, the exact accuracy, the analysed model and the PC. However, it can still serve as a rough estimation of the efficiency of the solvers. Even though Ghitti et al's solver exceed Wang et al's by an order of magnitude in run time, first_blood is even more efficient. It only requires 1:7 s to compute a heart cycle on average, roughly two times the real-time. The favourable computational cost creates the opportunity, for example, to produce a large database of virtual patients within days or to perform a comprehensive uncertainty analysis without the need for high-performance computers. Figure 9 shows how the relative error of physiological quantities (aortic diastole, aortic systole, cardiac output) converges to zero as the number of division points increases. Even with the coarser grid (only 1101 inner points overall), the relative errors are below 1%. By increasing the number of grid points, the error drops according to a power law function. The zeros of the y-axes are set as the limit of the fitted power law functions is zero. The run time increases according to a quadratic polynomial function (see bottom right of Figure 9). Although the number of grid points increases the computational time in space linearly, it also reduces the time step linearly. For a given simulation time, more calculated time steps are needed. This Subsection presents the results from the simulation of the human arterial circulation and the heart. Figure 10 depicts temporal history (the pressure and the flow rate) at five locations: aorta, carotid, radial, femoral, and anterior tibial artery. The figure also contains simulation data from first_blood with the presented model of the arteries and the heart, and data from the literature. 52 The amplitudes of the pressure, both the systolic and the diastolic ones, are in a physiologically acceptable range. Moreover, the signal shapes in various parts of the arterial system follow the trends. Regarding the velocity (or flow rate), the model contains an overshoot for the cardiac output, being 0.75 m/s (or 5.35 L/ min on average), while the literature suggests 0.4 m/s (or 4.57 L/min). However, the further pulse waves are correct both in shape and magnitude. Figure 11 shows the results for the heart, both the pressures and the flow rates. The vertical lines also indicate the state changes at the mitral and aortic valves, while the dashed blue curves are from the literature. 52 The pressure of the atrium remains constant since the 0D model of the heart prescribed it. As the ventricular pressure increases due to the change in the elastance, it reaches the aortic pressure, opening the aortic valve. The blood quickly starts to flow from the atrium to the ventricle. Once the excitation mitigates, the ventricular pressure drops below the aortic, closing the valve. A similar process can be observed with the mitral valve and the atrium pressure. Since ideal diodes model both heart valves, the model cannot capture the backflow effect. The dashed lines indicate literature data, 53 and there are some discrepancies.
Due to the simple left heart model, some phenomena cannot be modelled. The atrium pressure is constant, that is, it cannot follow real-life variation. The average flow rate of the heart slightly overshoots the literature data, but this may be due to the aforementioned lack of backflow. One possible solution is the introduction of hysteresis to the diodes F I G U R E 1 1 Results in the heart model, pressure ( p) above and flow rate (q) below as a function of time (t) for one cycle. The figure also indicates the state changes at the valves. Solid lines are simulation results from first_blood, and the dashed lines are from the literature: pressure data from Bucelli et al 53  modelling the heart valves. Again, the question of the correct choice of parameters to achieve the physiologically reasonable solution emerges. The direct measurement of this phenomenon is not feasible, and the literature lacks modelling of such features. Additional parameter calibration or model extensions are inevitable to increase the accuracy of the model further. Even though the model is low-dimensional with low computational time, it contains more than 1000 independent parameters; thus, the question is how to decrease their number. Jones et al 6 proposed a grouping based on the location of the arteries, then assigned the same multiplier factor to each parameter within the group. The question is whether this (or any other) grouping technique reduces the size of the search space, thus decreasing the accessibility of the optimum.

| SUMMARY
This paper introduced a 1D-0D combined, general solver called first_blood for hemodynamic phenomena. The main features are the following. first_blood is • computationally efficient, the simulation time of the whole arterial system with 99 vessels is only 2-3 times larger than the real-time, that is, the calculation of one heartbeat requires only 2 s. • The code is open source on GitHub. 47 • Lumped models can consist of linear (resistor, capacitor, etc.) and nonlinear 0D elements (e.g., valve, diode). The heart and the peripherals are handled with 0D models in this case study. • Finally, first_blood is modular in topology. The solver is general since the topology of the arteries, and the lumped models can be arbitrarily changed, and any additional organs or circulatory parts can be easily added. Moreover, the applicability of the solver demonstrated results corresponding to the arterial system, the heart, and the peripherals. The main arterial branches are resolved in 1D, and the governing fluid mechanical equations are solved using the method of characteristics. The material of the vessels is considered viscoelastic and is described with the Poynting-Thomson model.
The critical point is to set the model parameters correctly. They are determined partly based on the data from the literature and partly by manual tuning. The overall quality of the simulation data is acceptable, as it gives a good approximation of literature results, even without sophisticated calibration of the parameters. Moreover, the model extensions are available, as the combined 1D-0D solver is general, and any model topology can be accurately solved. Pulmonary, venous or coronary circulation elements can be added to further research to increase the accuracy of the model and the reliability in a physiological aspect.