The model contains a system of ordinary differential equations (ODEs):

Where X(t) and Y(t) are the system variables determined by the functions F : R^{n + m + l + 1} to Rⁿ and G : R^{n + m + l + 1} to R^m, respectively, and P denotes the system parameters. Note that Y can be expressed in terms of X and is selected into a separate set only for a clearer physiological interpretation of the model.

The system (1) describes continuous behavior of the model over time t and is accompanied by a number of discrete events corresponding to instantaneous changes in the model dynamics (for example, transition from systole to diastole). The event consists of a trigger and assignments. The trigger is the logical function T : R^n+m+l+1 → {true, false}. The event is considered triggered at the time point t’ if the value of this function changes from false to true at t = t′:

The event assignments are defined by the functions A_X : R^n+m+l+1 to Rⁿ and A_P : R^n+m+l+1 to R^l describing changes in X and P, respectively:

If any event is triggered at the time point t′, then solving of the system (1) is automatically stopped, the event assignments (2) are performed, and solving of the Cauchy problem is restarted from t = t′ with new initial values X(t′) = X′ and new parameter values P   = P′.

Such models, coupling the continuous and discrete approaches, are called the hybrid models (Stéphanou and Volpert, 2016). Parameters of the Cauchy problem (1) changed by the events (2) constitute the set of the model variables together with X and Y.

We consider the hybrid model to be in an equilibrium state at t = t_SS, if values of some variables Q₁, …, Q_k (called equilibrium variables) does not change in time: ∀t > t_SSQ_i(t) = Q_i(t_SS), i = 1, …, k. It is clear that not every such model with discrete events would have an equilibrium state. For instance, in chaotic triggering of some event, holding the equality F(X, Y, P, t) = 0 at t = t_SS does not guarantee its holding for all t > t_SS. In the case of modeling the cardiovascular system, discrete events represent switching between systole and diastole stages. Therefore, the model tends to fall into periodic behavior with a period equal to the cardiac cycle length. Thus, while most variables in X have non-zero derivatives, some variables in P changed by events are in equilibrium. For example, systemic arterial pressure is the dynamic variable dependent on differential equations. This variable increases in systole and decreases in diastole. At the same time, values of systolic and diastolic blood pressures are calculated at the moment of switching between those two stages. These quantities are in equilibrium (their value does not change in time).

To find numerical solutions of the direct problems, we applied the VODE solver (Brown et al., 1989; Cohen and Hindmarsh, 1996) supporting the automatic detection of time points at which discrete events are triggered.

Modularity could be considered as a principle of biological organization (Hartwell et al., 1999; Alon, 2003). Therefore, a modular approach to the modeling of complex biochemical systems has been actively developing in the last few years (Blinov et al., 2008; Hernández et al., 2009; Neal et al., 2014).

We define a module (Figure 1) as a part of a mathematical model describing a particular biological subsystem formulated as a separate block and integrated with the rest of the model using an explicit interface. Generally, modules can be treated as separate models with arbitrary mathematical formalism (ODE, stochastic, agent-based model, etc.) aggregated into modular models representing the larger biological systems. Such modules can also contain modular models, thus, forming a hierarchical structure with several nesting levels of modules. Every module defines interface variables used to connect modules with each other. In our case, the interface is defined using mathematical variables and parameters of the module. The input variables serve as the module parameters, these values must be calculated in other parts of the model and then passed to the module. The output variables are calculated inside the module but can be used outside it. Established connections in the modular model show which variables should be passed from one module to another. In the current study, we consider only modules containing sets of ODEs and discrete events (i.e., hybrid models described earlier). In that case, the modular model can be transformed into a “flat” hybrid model with the same formalism by aggregating all equations and events from all modules and resolving connections. For more details, see Kutumova et al. (2012) and Kiselev and Kolpakov (2013).

When modules function in significantly different time scales, an agent-based approach can be used to optimize the calculations. Applications of this approach span a broad range of areas from modeling the adaptive immune system to predicting the spread of epidemics (Macal and North, 2010). An agent is an autonomous entity which acts independently according to certain rules and interacts with other agents. Essentially, the agent is a black box that receives a signal, sends a response and has specific mathematical formalism and numerical methods inside.

When the agent-based approach is used for simulation of the modular models, each module is considered as the separate agent consisting of some mathematical model, numerical solver and time span (initial time, time step, and completion time). Coordination between agents is provided by a scheduler. The agent step corresponds to the model simulation from the current span point to the next one. Before a new step, the scheduler notifies the agent about last changes made to its input interface variables. When the step is finished, the agent sends to the scheduler the changes that it made to the output interface variables. Note that the agent spans are used only to determine time points of the agent interactions, while numerical calculations inside each agent involve automatic computation of time steps necessary for required accuracy.

The agent-based simulation goes as follows (scheme 1).

The problem arises when two agents have drastically different time scales. For example, we have the “Cardiovascular system” module describing very fast processes of the heart pumping and blood flow across the vascular system, and the “Renal system” module with the long-term regulation of salt-water balance, total blood volume and hormone levels. If we combine both modules into one ODE system, we get a very stiff problem. When using the agent-based form of the model, we still have the problem related to the simulation of the “fast” model over a very large time span, which is very time consuming. However, if the “fast” model has an equilibrium state, we can use another approach to simulation of the agent-based models.

Let’s consider the model with two agents shown in Figure 1. For numerical simulation we will use time span t₀, t₁, t₂, …, t_N, where t_{i + 1} = t_i + Δ, i = 1, …, N. At each time point t_i, the agent corresponding to Module 1 calculates a new value y(t_i) and sends it to Agent 2 where it is used as a new value for the parameter a. Then, Agent 2 starts the next simulation step with the value a = y(t_i). Exact solution of the agent-based model would be the exact solution of the ODE system represented by the “flat” version of the model. Thus, using the agent based approach is equivalent to solving both equations separately between the agent exchanges. That means that when solving the equation for x at each time step, we use the numerical value of y from the previous time point and introduce the local error O(Δ²). Let x¯(t) and y¯(t) be the simulation results of the flat model, whereas x(t) and y(t) are the simulation results of the agent-based model. Suppose, we have no error at the time point t₀:

Using the Taylor series for both solutions, we obtain:

When considering that y¯(t0)=a, we derive:

Now let’s suppose that T_x≫T_y and Agent 2 has an equilibrium point for all values of parameter a, i.e., if at some t = t_i, the model is in equilibrium x(t_i) with a = a₁, then changing the value of a to a₂ and solving the Cauchy problem for x with the new initial value equal to x(t_i) gives new equilibrium x(t_i + t′). An example result of such model simulation is presented in Figure 2. One can see that after each exchange between agents, the “fast” agent rapidly approaches equilibrium. Therefore, this agent simulation until the next moment of exchange becomes excessive and unnecessary.

A way to optimize scheme 1 for the numerical calculation in this particular case is to stop simulation for the “fast” agent after it reaches the equilibrium at each agent step. Thus, at each time point t_i, we set:

Where t′ is the time interval during which Agent 2 reaches the equilibrium after the exchange at the time point t_i–1.

There is also a way to get a smaller error and simultaneously keep a large time step between agent exchanges. Obviously, the time step of agents Δ is much larger than the time interval t′. Therefore, we would introduce a smaller error if we use a new equilibrium value for the previous span point rather than the next one. Instead of (3), we set:

Applying both described optimizations, we get scheme 2 for the agent modeling:

Comparison between both schemes of the simulation is given in Figure 3.

Parameter estimates performed in this work for the model calibration were based on minimization of the distance function defined as a normalized sum of squared differences (Hoops et al., 2006) between simulated equilibrium values Q₁, …, Q_k and clinical measurements Q1exp,…,Qkexpof physiological quantities:

Where ω_min = min_iω_i and weights ω_i are used to make all quantities to have similar importance.

To keep certain model variables W₁(t), …, W_p(t) in physiological limits, we considered known constraints Wimin≤Wi(t)≤Wimax and additionally minimized the penalty function (Runarsson and Yao, 2000):

Calculated for t ∈ [t_ss, t_ss + d], where t_SS is the model equilibrium time point and d is the length of the cardiac cycle. Note that, as mentioned above, though the model is in equilibrium, some of its variables (e.g., arterial pressure) are dynamic and change their values during the cardiac cycle due to the presence of discrete events.

The process of parameter fitting was based on the stochastic ranking method suitable for the constrained optimization (Runarsson and Yao, 2000).

BioUML (homepage)² is an integrated Java platform for modeling of biological systems (Kolpakov et al., 2019). It supports a comprehensive range of tools for systems biology, including visual modeling, simulation, parameter estimation and a number of numerical methods. Key features of the software used in this work:

A visual approach to mathematical modeling involves the creation and work with mathematical models as diagrams. Thus, each element of the model (equation, event, interface variable, agent, connection) corresponds to an element of the diagram: edge or node. The visual modeling implies using some kind of formal graphical notation, so the visual representation of each particular element depends on their mathematical properties. A common standard for the graphical notation in systems biology is SBGN (Systems Biology Graphical Notation, Le Novère et al., 2009). However, it is mostly used for mathematical models comprising pathways of processes (e.g., biochemical reactions) and it lacks the elements representing the arbitrary differential or algebraic equations. Thus, in the current study, we use the visual notation developed in the BioUML platform. The visual representation makes the inner structure of the model more explicit and facilitates understanding of the model.

The comprehensive model of the human cardiovascular and renal systems is based on the range of models of renal hemodynamics (Karaaslan et al., 2005, 2014; Hallow et al., 2014; Hallow and Gebremichael, 2017a) originating from the Guyton model (Guyton et al., 1972), and the model of blood circulation system (Proshin and Solodyannikov, 2006). Below we provide a brief description of the modular representation, while details of the equations, parameters and variables for each module are given in the Supplementary Appendix Tables A1–A4. Totally, the model contains 20 modules, 10 discrete events, 185 equations, 132 parameters, and 160 variables.

On the top level, the model can be introduced as an interaction of two main agents: the renal dynamics is determined in minutes, whereas the cardiovascular processes take a fraction of a second (Figure 4). These agents have five main connections corresponding to physiological quantities that are calculated in one agent and directly affect the dynamics of another one. In the cardiovascular system, such quantities are mean arterial pressure, cardiac output and a value of hematocrit. The renal system is responsible for regulation of the total body blood volume and concentration of angiotensin II bound to the AT1 receptors.

For clarity, the model can be further decomposed into 20 functional modules: 11 in the cardiovascular sub-model and 9 in the renal sub-model (Figure 5). Six main cardiovascular modules form a circular system of compartments (Proshin and Solodyannikov, 2006): Left ventricle, LV (denoted HL in the model notation), Systemic arteries (AL), Systemic veins (VL), Right ventricle, RV (HR), Pulmonary arteries (AR), and Pulmonary veins (VR). Each i-th compartment, i ∈ {HL, AL, VL, HR, AR, VR}, is characterized by pressure P_i, volume V_i, unstressed volume ω_i, and elastance G_i. Blood flow F_ij between the i-th and j-th compartments is determined by the difference of their pressures, and can be simply written as:

Where Y_ij is conductivity of the respective part of the circulatory system. Thus, an important variable affecting blood flow from arteries to veins is conductivity of microvessels (arterioles, capillaries and venules). We allocated the calculation of this variable for systemic and pulmonary circulation into separate modules “Systemic microvessels” and “Pulmonary microvessels,” respectively.

One of the main ideas of the model constructed by Proshin and Solodyannikov (2006) is the use of discrete events that determine the instantaneous change in the parameters of the left and right ventricles at the moment of a cardiac cycle start (diastole to systole transition), and at the moment when blood ejection from the ventricles stops. Since the LV and RV ejection time is different (Leighton et al., 1971; Hirschfeld et al., 1975), the model contains 2 discrete events of systole to diastole transition in the LV and RV modules, respectively. Each of these 2 events, in particular, changes the value of the corresponding indicator Systole_L or Systole_R from 1 (ejection is in progress) to 0 (ejection is finished). The common indicator Systole tracing the moment when blood ejection from both ventricles comes to the end is placed into the module “Coordination of systolic time intervals.” It is equal to 0 if Systole_L = 0 and Systole_R = 0, and 1 otherwise.

A model of oxygen exchange describing oxygen consumption in tissues and oxygen debt based on the amount F_ALVL was taken unchanged from the work (Proshin and Solodyannikov, 2006) and inserted into the module “Tissue metabolism.”

The last module “Neurohumoral control” in the heart sub-model defines the cardiac center as a control system forming the output value of a hypothetical neurohumoral factor H from the sum signal of receptor activities. The model by Proshin and Solodyannikov (2006) considers four types of receptors, namely stress (nS), weariness (nD), and respiratory receptors (nV), as well as baroreceptors (nB). Receptors respond to various internal factors of the body and external influences and transmit to the cardiac center the signals calculated by the formulas:

Where P_AL is systemic arterial pressure (determined in the module “systemic arteries”), VO₂ and DO₂, respectively denote venous oxygen content and oxygen debt (“tissue metabolism”), st ∈ [0,1] define the level of steroid hormones in the blood, whereas r⁻(φ) = 1−r⁺(φ) and r⁺(φ) is the sigmoid function:

The initial values of the constants α, β, and φ₀ for all receptors are listed in the study (Proshin and Solodyannikov, 2006).

The renal sub-model incorporates the renin-angiotensin-aldosterone system (RAAS) pathway divided into three main modules:

Dynamics of these modules depends on the renal sympathetic nerve activity calculated in the module “Nervous system,” whereas concentration of AT1-bound angiotensin II directly effects on tubular sodium reabsorption (“Sodium”) and renal vascular resistance (composed of resistances of afferent/efferent arterioles, and interlobar/arcuate/interlobular arteries), as well as renal blood flow and filtration (“Glomerular filtration”). The remaining modules are responsible for the calculation of atrial natriuretic peptide and antidiuretic hormone concentrations (“Hormonal system”), urine flow (“Diuresis”), and volumes of total body water, blood, and extracellular fluid (“Body fluids”).

Implementation of the model based on parts of the existing models required us to make a number of changes to them with the following goals:

Below we provide a description of the primary physiological processes involved in the model, while the detailed formulas can be found in Supplementary Appendices A (Table A2), B.

Each equilibrium parameterization of the model can be considered as a unique virtual patient (Cheng et al., 2017). To deal with a variety of such patients (which is related to the search for different equilibriums of the model depending on the starting values), we need to have a tool for scaling of physiological variables. For this purpose, Proshin and Solodyannikov (2006) used the general scaling scheme representing the dependence of biological variables a_i (elasticity of vascular walls, conductivity of vessels, body oxygen demand, etc.) on body mass m by an allometric power-laws (West et al., 1997; Proshin and Solodyannikov, 2006):

Where b_i are the scaling exponents and c_i are the normalization constants. In addition to this approach, we took into account the approximate distribution of blood (in percentage of total blood V) in the different parts of the circulatory system: 84% – systemic circulation (including 13% in systemic arteries) (Hall, 2011), 10% – pulmonary circulation (Gazioglu and Yu, 1967), 6% – heart. Pulmonary vessels include arterial (35% of pulmonary circulation), venous (45%) and capillary (20%) volumes (Gazioglu and Yu, 1967). Since in the model by Proshin and Solodyannikov (2006) pulmonary veins and capillaries are included in one compartment of the circulatory system, assuming equality of the starting V_HL and V_HR, we get the following formulas for calculating the initial volumes of the compartments:

The corresponding unstressed volumes we determined as:

Where estimated constants k_AL, k_VL, k_AR, k_VR ∈ [0.7, 1.0], and k_HL, k_HR ∈ [0.0, 0.3] (see the Supplementary Appendix Table A3 for details of such intervals selection).

Many parameters of the model correspond to physiological quantities, which values can be evaluated by the laboratory measurements. Thus, the model calibration consists in finding an equilibrium state satisfying the set of physiological constraints. The explicit form of these constraints depends on the studied problem and the status of the simulated patient (healthy or sick). For example, the normal range of the systolic/diastolic pulmonary artery pressure is 15–30/4–12 mmHg (Pagani et al., 1988). In patients with class II heart failure according to the New York Heart Association classification, the corresponding values are higher: 35.4 ± 8.8/14.4 ± 5.8 mmHg (Murch et al., 2015), whereas in patients with pulmonary arterial hypertension these values can achieve 84.0 ± 23.0/37.0 ± 13.0 mmHg (Gan et al., 2006).

In this work, we present the model’s ability to simulate healthy subjects as well as patients with the most common cardiovascular diseases. We considered the following test cases:

The model calibration included two main steps.

The resulting equilibrium values of the model within the ranges of a healthy person are listed in the Supplementary Appendix Tables A3, A4. These values are also provided in the BioUML web-implementation of the model. Physiological quantities of the equilibrium states representing abstract patients 1–5 are given in Table 2 and introduced as the model states stored as separate documents and applicable to the model (see the Availability section below).

As examples of using the model, we consider two standard problems faced by researchers in the study of human physiology: comparison of different physiological states and comparison of sensitivity of different patients to the change in physiological conditions. To produce the comparative plots, we used the capabilities of the Jupyter notebook embedded in BioUML. This application is designed to create and share documents that contain live code, equations, visualizations and narrative text. The Jupyter files comprising the implementation of both use cases described below are available in the web-version of our software (see the Availability section).