Thoroughly Calibrated Modular Agent-Based Model of the Human Cardiovascular and Renal Systems for Blood Pressure Regulation in Health and Disease

Here we present a modular agent-based mathematical model of the human cardiovascular and renal systems. It integrates the previous models primarily developed by A. C. Guyton, F. Karaaslan, K. M. Hallow, and Y. V. Solodyannikov. We performed the model calibration to find an equilibrium state within the normal vital sign ranges for a healthy adult. We verified the model’s abilities to reproduce equilibrium states with abnormal physiological values related to different combinations of cardiovascular diseases (such as systemic hypertension, chronic heart failure, pulmonary hypertension, etc.). For the model creation and validation, we involved over 200 scientific studies covering known models of the human cardiovascular and renal functions, biosimulation platforms, and clinical measurements of physiological quantities in normal and pathological conditions. We compiled detailed documentation describing all equations, parameters and variables of the model with justification of all formulas and values. The model is implemented in BioUML and available in the web-version of the software.

Introduction

Mathematical models provide tools for understanding of human physiology via integration and analysis of biological data from multiple ranges and time scales under normal and pathological conditions. A number of such models have already been used to investigate different aspects of individual systems and processes of the human body (Guyton et al., 1972; Ikeda et al., 1979; Uttamsingh et al., 1985; Ottesen et al., 2004; Karaaslan et al., 2005, 2014; Proshin and Solodyannikov, 2006; Abram et al., 2007; Thomas et al., 2007; Hester et al., 2011; Paeme et al., 2011; Hallow et al., 2014, 2021; Hallow and Gebremichael, 2017a; Rosalina et al., 2019). Most of them originate from the control-theory block model of circulatory regulation proposed by Guyton et al. (1972). This model was implemented in different programming languages (Moss et al., 2012) and reused in further whole-body models (Abram et al., 2007; Hester et al., 2011). Guyton et al. (1972) established the foundation for our understanding of the relations between long-term blood pressure level and sodium and water balance. And in this sense, their model is innovative and revolutionary (Montani and Vliet, 2009; Beard, 2018). However, its application for research purposes is difficult due to the model complexity, the lack of a complete detailed description, physiological limitations for some parameters, and some misconceptions (Montani and Vliet, 2009).

To solve narrowly focused problems, more compact models have been created recently. Thus, Karaaslan et al. (2005) explored the renal sympathetic nerve activity related mechanisms affecting kidney functions and causing increase of arterial pressure in hypertension. Hallow et al. (2014) simulated pharmacodynamic effects of the antihypertensive therapy, and then created a detailed model of renal physiology (Hallow and Gebremichael, 2017a) applied for investigation of salt-sensitive hypertension (Hallow and Gebremichael, 2017b). A number of models focused on cardiovascular hemodynamics (Proshin and Solodyannikov, 2006; Thomas et al., 2007; Paeme et al., 2011; Rosalina et al., 2019).

Accumulation of the particular models of human physiology leads to the next logical step: “gluing” of these models (Kassab and Guccione, 2019). However, this step results in the creation of complex entities. The only way to deal with such complexity is dividing the biological systems into subsystems, and the corresponding models into modules. This approach simplifies the modeling process and, at the same time, allows you to gradually improve the model by replacement and addition of modules. Another problem that arises in the model integration is related to the combination of processes acting at significantly different time scales. For example, a single heartbeat takes less than a second, whereas a renal function is associated with the long-term regulation of the salt and water balance. Merging the equations for these processes gives a very stiff problem which in the specified case can be solved by application of an agent-based approach.

At the present time, there are several biosimulation platforms designed for modeling and analysis of the human physiological processes. For instance, the Entelos PhysioLab platforms are implemented for drug development and have been applied to projects in diabetes, rheumatoid arthritis, asthma, and skin sensitization (Rullmann et al., 2005; Maxwell and Mackay, 2008; Klinke, 2015). These projects are widely used by pharmaceutical companies (Pfizer, Novartis, etc.) in preclinical trials of new drugs, and capture the modeled physiology using a modular approach, assuming that smaller scale models are defined based on isolated components (cells or signaling pathways) connected together to reflect higher (intracellular) level behavior (Klinke, 2015). Another major project is the open-source, full-body human physiology engine BioGears (McDaniel et al., 2019). Its purpose is to provide realistic and comprehensive simulations for medical research and education. BioGears may be used as a standalone application or integrated with simulators, sensor interfaces, and other existing models. According to information from the developers, this platform operates with lumped parameters and, therefore, is not aimed at creating personalized models, but uses “average” person for analysis. One more environment for modeling and simulation of integrative human physiology is the HumMod software (Hester et al., 2011). It describes cardiovascular, respiratory, renal, neural, endocrine, skeletal muscle, and metabolic systems, and is constructed from empirical data obtained by authors from scientific literature. As in the case of BioGears, the user of HumMod defines a number of patient basic parameters (gender, height, weight, etc.) and operates with an average patient in the normal or pathological state. As an example of the environment for modeling of the circulatory system, we want to single out the Samara-Dialog platform designed on the basis of the cardiovascular hemodynamics model proposed by Proshin and Solodyannikov (2006)¹. This environment is primarily intended to study the status of an athlete in the training process, but can also be used to simulate a wide range of the human cardiovascular system pathologies, including heart arrhythmias, ventricular dysfunctions, valvular failure, hypertension, etc.

All of the mentioned models and simulation platforms were created to investigate mechanisms of different cardiovascular and renal diseases in several abstract conditions, where parameter values are fixed on the basis of some average normal or pathological evaluations. However, cardiac models should account for the fact that humans vary (Wiśniowska et al., 2017). If we consider, in addition, the mathematical modeling to predict real-world effectiveness of drug interventions, we find that the most important limitation is the lack of external validation applying other data than those used for developing the models. In their review, Panayidou et al. (2016) concluded that such modeling is not widely used at present and not well validated.

Thus, in this work, we focus not only on creating a modular agent-based model of the cardiovascular and renal systems, but also perform detailed validation of parameters in order to conform to physiological ranges and reproduce equilibrium states corresponding to various combinations of cardiovascular diseases in real patients. This can later be used for personalized modeling and individual predictions on a case-by-case basis. For operation with the model, we used the BioUML software, which provided all the necessary tools for our research, and which our team has been developing since 2002 (Kolpakov et al., 2019).

Materials and Methods

Mathematical Formulation

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

$$\begin{array}{c}\{\begin{array}{c}\frac{\mathit{dX}}{\mathit{dt}}=\mathit{\text{F}}(\mathit{X},\mathit{Y},\mathit{P},\mathit{t}),\hfill \\ \mathit{\text{X}}(\mathit{0})={\mathit{\text{X}}}_{\mathbf{0}},\hfill \\ \mathit{\text{Y}}(\mathit{t})=\mathit{\text{G}}(\mathit{X},\mathit{Y},\mathit{P},\mathit{t})\hfill \end{array}\hfill \\ \mathit{\text{X}}=\left(\begin{array}{c}\hfill {\mathit{\text{X}}}_{\mathit{1}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\mathit{\text{X}}}_{\mathit{n}}\hfill \end{array}\right),\mathit{Y}=\left(\begin{array}{c}\hfill {\mathit{\text{Y}}}_{\mathit{1}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\mathit{\text{Y}}}_{\mathit{m}}\hfill \end{array}\right),\mathit{P}=\left(\begin{array}{c}\hfill {\mathit{\text{P}}}_{\mathit{1}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\mathit{\text{P}}}_{\mathit{l}}\hfill \end{array}\right)\hfill \end{array}(1)$$

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′:

$$\begin{array}{c}\exists \mathrm{\delta }>\mathbf{0}:\forall \mathit{\text{t}}\in [{\mathit{\text{t}}}^{{}^{\prime }}-\mathrm{\delta },{\mathit{\text{t}}}^{{}^{\prime }})\mathit{\text{T}}(\mathit{\text{X}}(\mathit{t}),\mathit{\text{Y}}(\mathit{t}),\mathit{P},\mathit{t})\hfill \\ =\mathit{false}\&\mathit{\text{T}}(\mathit{\text{X}}({\mathit{\text{t}}}^{{}^{\prime }}),\mathit{\text{Y}}({\mathit{\text{t}}}^{{}^{\prime }}),\mathit{P},{\mathit{\text{t}}}^{{}^{\prime }})=\mathit{\text{true}}\hfill \end{array}$$

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:

$$\{\begin{array}{c}{\mathit{\text{X}}}^{\prime }={\mathit{\text{A}}}_{\mathit{X}}(\mathit{\text{X}}({\mathit{\text{t}}}^{\prime }),\mathit{\text{Y}}({\mathit{\text{t}}}^{\prime }),\mathit{P},{\mathit{\text{t}}}^{\prime }),\hfill \\ {\mathit{\text{P}}}^{\prime }={\mathit{\text{A}}}_{\mathit{P}}(\mathit{\text{X}}({\mathit{\text{t}}}^{\prime }),\mathit{\text{Y}}({\mathit{\text{t}}}^{\prime }),\mathit{P},{\mathit{\text{t}}}^{\prime })\hfill \end{array}(2)$$

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.

Modular Modeling

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).

FIGURE 1

Figure 1. An example of the modular representation of a two-equation model. Module 1 includes one ODE for y (output interface variable), as well as two parameters b and T_y. Module 2 comprises one ODE for x, one input interface variable a, and one parameter T_x.

Agent-Based Modeling

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).

(1) Pick the agent with the least current model time value.

(2) Pass to the selected agent the changes done by other agents to its input interface variables.

(3) Perform the step in time for the selected agent (i.e., run simulation of the mathematical model incorporated in the agent).

(4) Send the changes made by the selected agent to its output interface variables to the scheduler (for further translation to other agents).

(5) Check if the simulation is finished. If not, go to 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 $\overline{x}\left(t\right)$ and $\overline{y}\left(t\right)$ 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₀:

$$\mathit{\text{x}}({\mathit{\text{t}}}_{\mathbf{0}})=\overline{\mathit{\text{x}}}({\mathit{\text{t}}}_{\mathbf{0}}),\mathit{\text{y}}({\mathit{\text{t}}}_{\mathbf{0}})=\overline{\mathit{\text{y}}}({\mathit{\text{t}}}_{\mathbf{0}})=\mathit{\text{a}}$$

Using the Taylor series for both solutions, we obtain:

$$\overline{\mathit{\text{x}}}({\mathit{\text{t}}}_{\mathbf{0}}+\mathrm{\Delta })=\overline{\mathit{\text{x}}}({\mathit{\text{t}}}_{\mathbf{0}})+\mathrm{\Delta }{\mathit{\text{T}}}_{\mathit{X}}(\overline{\mathit{\text{y}}}({\mathit{\text{t}}}_{\mathbf{0}})-\overline{\mathit{\text{x}}}({\mathit{\text{t}}}_0))+\mathit{\text{O}}({\mathrm{\Delta }}^{\mathbf{2}})$$

$$\mathit{\text{x}}({\mathit{\text{t}}}_0+\mathrm{\Delta })=\mathit{\text{x}}({\mathit{\text{t}}}_0)+\mathrm{\Delta }{\mathit{\text{T}}}_{\mathit{X}}(\mathit{\text{a}}-\mathit{\text{x}}({\mathit{\text{t}}}_{\mathbf{0}}))+\mathbf{\text{O}}({\mathrm{\Delta }}^{\mathbf{2}})$$

When considering that $\overline{y}\left(t_0\right)=a$, we derive:

$$\mathit{\text{x}}({\mathit{\text{t}}}_{\mathbf{0}}+\mathrm{\Delta })=\overline{\mathit{\text{x}}}({\mathit{\text{t}}}_0+\mathrm{\Delta })+\mathit{\text{O}}({\mathrm{\Delta }}^{\mathbf{2}})$$

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.

FIGURE 2

Figure 2. A simulation example for the model defined in Figure 1 with T_x≫T_y and time Δ = 1000 between agent interactions.

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:

$$\mathit{\text{x}}\left({\mathit{\text{t}}}_{\mathit{i}}\right)=\mathit{\text{x}}\left({\mathit{\text{t}}}_{\mathit{i}-\mathbf{1}}+{\mathit{\text{t}}}^{{}^{\prime }}\right)(3)$$

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:

$$\mathit{\text{x}}\left({\mathit{\text{t}}}_{\mathit{i}}\right)=\mathit{\text{x}}\left({\mathit{\text{t}}}_{\mathit{i}}+{\mathit{\text{t}}}^{{}^{\prime }}\right)$$

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

(1) Pick the agent with the least current model time value. “Fast” agents have a lower priority and, thus, should perform their steps after “slow” agents.

(2) If the selected agent is “slow,” use scheme 1.

(3) Pass to the selected agent the changes done by other agents to its interface variables (For the model in Figure 1: a   = y_i).

(4) Perform the step in time for the selected agent (i.e., run simulation of the mathematical model incorporated in the agent).

(a) If the agent is “slow”, perform the numerical calculations until the model reaches the next time step x(t_i) → x(t_{i + 1}) = x(t_i + Δ).

(b) If the agent is “fast”, perform the numerical calculations until the model reaches the equilibrium x(t_i) → x(t_i + t′). Set new value for the current time point x(t_i) = x(t_i + t′). Set the model time of the agent to the next span point t = t_{i + 1} (thus, the value of x(t_{i + 1}) will be calculated at the next step).

(5) Send changes made by the selected agent to its interface variables to the scheduler (for further translation to other agents).

(6) Check if the simulation is finished. If not, go to 1.

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

FIGURE 3

Figure 3. Comparison between two schemes of the two-agent model simulation. The second scheme, instead of the long and unnecessary numerical calculations of the “fast” Agent 2, searches for the equilibrium at each point of the time span.

Parameter Estimation

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 $Q_1^{exp},\mathrm{\dots },Q_k^{exp}$of physiological quantities:

$${\mathit{\text{f}}}_{\mathit{dist}}=\sum _{\mathit{i}=\mathbf{1}}^{\mathit{k}}\frac{{\mathrm{\omega }}_{\mathit{min}}}{{\mathrm{\omega }}_{\mathit{i}}}{\left({\mathit{\text{Q}}}_{\mathit{i}}\left({\mathit{\text{t}}}_{\mathit{ss}}\right)-\mathit{\text{Q}}_{\mathit{i}}{}^{\mathit{exp}}\right)}^2,{\mathrm{\omega }}_{\mathit{i}}={\mathit{\text{Q}}}_{\mathit{i}}^{\mathit{exp}}(4)$$

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 $W_i^{min}\le W_i\left(t\right)\le W_i^{max}$ and additionally minimized the penalty function (Runarsson and Yao, 2000):

$$\begin{array}{c}{f}_{\mathrm{penalty}}={\displaystyle \sum _t({\displaystyle \sum _{i=\mathit{1}}^P\max \{\mathbf{0},W_i^{\min }-W_i(t){\}}^{\mathbf{2}}}}\hfill \\ +{\displaystyle \sum _{i=\mathbf{1}}^p\max \{\mathbf{0},W_i(t)-W_i^{\max }{\}}^{\mathbf{2}}})\hfill \end{array}(5)$$

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).

Modeling Platform

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:

• The opportunity to work independently in the local (standalone) version of the program or through the web interface in collaboration with other researchers.

• Plugin-based architecture of the platform allowing to design new types of models and to implement required methods for the numerical analysis of them.

• An editor for visual modeling of biological systems using modular and agent-based approaches.

• The embedded VODE solver (Cohen and Hindmarsh, 1996) ported to Java and suitable for the numerical simulation of hybrid models with ODE systems and discrete events.

• A number of embedded methods (in particular, the stochastic ranking evolution strategy (Runarsson and Yao, 2000) which we preferably used in this study) for the model parameter fitting based on the reference ranges of physiological parameters in normal and pathological states.

• Integration with the JupyterHub³ for interactive data analysis.

• Support of the SBML standard (Hucka et al., 2019) for model exchange.

Visual Modeling

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.

Results

Modular Structure 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.

FIGURE 4

Figure 4. Modular modeling of the cardiovascular and renal systems. The original model decomposed into two modules in accordance with human physiology. The “Cardiovascular system” module is responsible for calculation of such variables as mean arterial pressure (MAP) and cardiac output (CO). The “Renal system” module in turn regulates total body blood volume (V) and concentration of angiotensin II bound to the AT1 receptors (AT1-ANG II). The value of hematocrit (Hct) is constant. It is determined in “Cardiovascular system” and passed to “Renal system” for use.

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:

$${\mathit{\text{F}}}_{\mathit{ij}}={\mathit{\text{Y}}}_{\mathit{ij}}\left({\mathit{\text{P}}}_{\mathit{i}}-{\mathit{\text{P}}}_{\mathit{j}}\right)$$

FIGURE 5

Figure 5. Modular agent-based physiological model of the cardiovascular and renal systems. (A) An example of the “Pulmonary arteries” module representation. (B) The model implemented in the BioUML platform is divided into nine modules responsible for the kidney function (green) and 11 modules that simulate dynamics of the cardiovascular system (purple). Blue arrows indicate directed connections between the renal and cardiovascular sub-diagrams (Figure 4). For visual simplicity of the diagram, we added transitional nodes (busses) which are used for connections between modules. The busses corresponding to one variable can be located far apart in the diagram. A complete list of equations for each module is given in the Supplementary Appendix Table A2.

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:

$$\begin{array}{c}\mathit{nB}={\mathit{\text{r}}}^{-}({\mathit{\text{P}}}_{AL}),\mathit{nV}={\mathit{\text{r}}}^{-}\left({\mathit{\text{VO}}}_{\mathbf{2}}\right),\hfill \\ \mathit{nD}={\mathit{\text{r}}}^{+}\left({\mathit{\text{DO}}}_2\right),\mathit{nS}={\mathit{\text{r}}}^{+}\left(\mathit{st}\right)\hfill \end{array}$$

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:

$${\mathit{\text{r}}}^{+}(\mathrm{\varphi })={\mathit{\text{r}}}^{+}(\mathbf{\alpha },\mathbf{\beta },\mathbf{\varphi },{\mathbf{\varphi }}_{\mathbf{0}})=\frac{\mathbf{1}-\mathit{\text{exp}}\left(-\mathbf{\alpha }\left(\mathbf{\varphi }-{\mathbf{\varphi }}_{\mathbf{0}}\right)\right)}{\mathbf{1}+\mathbf{\beta }\cdot \mathit{\text{exp}}\left(-\mathbf{\alpha }\left(\mathbf{\varphi }-{\mathbf{\varphi }}_{\mathbf{0}}\right)\right)}$$

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:

• Secretion of renin (the “Renin” module).

• Generation of angiotensin I with consistent formation of angiotensin II, -IV, -(1–7), and activation of AT1/AT2 receptors (“Angiotensin”).

• Secretion of aldosterone (“Aldosterone”).

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”).

Details and Updates to the Base Models

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:

• Relation of these parts with each other.

• Inclusion of the clinically measurable variables to physiological ranges.

• Getting the ability to reproduce pathologies.

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.

Targets of Angiotensin II

Angiotensin II exerts physiologic actions via binding to receptors on cells of different organs (the kidneys, heart, blood vessels, etc.). The leading role in such actions belongs to the AT1 receptors. Based on their localization (Allen et al., 2000), we have the following angiotensin II targets in the model.

• Vascular Smooth Muscle Angiotensin II exerts vasoconstrictor effect (Hughes, 1998). Thus, if the concentration of AT1-bound angiotensin II (AT1_ANGII) increases, then resistances of afferent/efferent arterioles (R_aa and R_ea) and interlobar/arcuate/interlobular arteries (R_preglom) rise resulting in the renal blood flow decline. Functions of AT1_ANGII effects on R_aa, R_ea, and R_preglom were suggested by Hallow et al. (2014). We assumed that the concentration of AT1-bound angiotensin II in the systemic arterioles is the same as that in the kidneys.

• Cardiac Muscle Angiotensin II-induced positive chronotropic and slight positive inotropic effects were demonstrated in the isolated dog heart (Kobayashi et al., 1978). At the same time, the action of angiotensin II in humans is associated either with no change in heart rate or with a reduction that is much smaller than that produced by other vasoconstrictors (Reid, 1996). Thus, when modeling the response to antihypertensive therapies, the influence of angiotensin II on cardiac muscle must be taken into account to simulate decrease in blood pressure without increasing heart rate, which, for instance, is the standard effect of such RAAS inhibitors as aliskiren, losartan, and enalapril (Konstam et al., 1993; Kamishirado et al., 1997; González-Abraldes et al., 2001; Stanton et al., 2003; Parrinello et al., 2009; Natarajan et al., 2016). In this regard, we considered in the model two targets of angiotensin II: baroreceptors and stress receptors. Such a decision was based on the facts that angiotensin II resets the baroreflex control of heart rate to a higher pressure (Reid, 1996) and increases the release of norepinephrine from the atria (Brasch et al., 1993).

• Adrenal Glands Angiotensin II stimulates aldosterone synthesis in adrenal zona glomerulosa cells via binding to the AT1 receptors (Bandulik et al., 2015). A function for modeling the AT1_ANGII effect on the aldosterone secretion rate was taken from the models by Karaaslan et al. (2005) and Hallow et al. (2014).

• Kidney, Proximal Tubules An increase in the level of AT1_ANGII leads to the rise of fractional proximal sodium reabsorption. The corresponding formula was found in the model by Karaaslan et al. (2005).

• Kidney, Juxtaglomerular Cells In vitro studies have suggested the expression of the AT1 receptors by juxtaglomerular cells. In this regard, there is the concept that angiotensin II directly inhibits renin secretion through a negative feedback. Although this concept has not yet been systematically proven in vivo and can be doubted (Neubauer et al., 2018), the factors affecting renin secretion (macula densa sodium sensing and renal sympathetic nerve activity) are not be sufficient to account for the rise in renin that occurs with therapeutic blockade of the RAAS (Hallow et al., 2014). Thus, following the model by Hallow et al. (2014) we used the formula of the direct negative feedback of AT1-bound angiotensin II on renin secretion introduced in the study (Hallow and Gebremichael, 2017a).

• Kidney, Mesangial Cells Angiotensin II provides a reduction in the glomerular filtration coefficient K_FG (Blantz et al., 1976) via decrease in total filtering surface area because of mesangial cell contraction (Schmitt et al., 1998). Angiotensin II receptors on mesangial cells belong to the AT1 subtype (Ardaillou et al., 1999). Thus, for the calculation of K_FG, we considered a product of its normal value and a linear function expressing inversely proportional relationship between K_FG and the normalized concentration of AT1_ANGII.

Peak Flow Rates Through the Heart Valves

The rate of the LV filling through the mitral valve in diastole is characterized by two peaks (Boogers et al., 2011; Caudron et al., 2011; Zhang et al., 2019), which can be estimated from the time derivative of the LV volume by cardiovascular magnetic resonance (Maceira et al., 2006a). In normal subjects, the LV inflow is greatest immediately after opening of the mitral valve (early peak), while the left atrial contraction is responsible for smaller inflow (active peak) (Caudron et al., 2011). The moment of the mitral valve opening in the model is the LV transition to diastole, when the pressure in pulmonary veins becomes greater than the pressure in the ventricle. The left atrial contraction can be associated with a positive value of the left atrial pulse wave, proposed in the model (Proshin and Solodyannikov, 2006). The right ventricular filling through the tricuspid valve has similar dynamics (Maceira et al., 2006b). As for the transaortic and transpulmonary flows, they reach their maximum values at the transition from diastole to systole. Therefore, calculation of the peak transvalvular flow rates can be introduced by the formalism of discrete events.

Physiological Quantities

In our work, we strived to create a model focused not on average values of physiological parameters, but on a variety of values including reference intervals for healthy people and possible pathological deviations from the norm. Since the baseline models did not accept these ranges as valid for some variables, we used experimental formulas obtained in population studies for the following variables.

• Plasma sodium concentration is calculated by the formula derived in the study (Nguyen and Kurtz, 2003) and obtained on the basis of the Edelman equation (Edelman et al., 1958).

• Antidiuretic hormone concentration is defined depending on serum osmolality (Hammer et al., 1980), which value can be considered as function of plasma sodium concentration and levels of glucose and urea in blood (Dorwart and Chalmers, 1975; Bhagat et al., 1984).

• Body fluids: For the calculation of the total blood volume and the extracellular fluid volume we applied the regression functions on the total blood water (Moore, 1967).

• Tubular water reabsorption rate depends on the glomerular filtration rate. In the absence of significant amounts of poorly reabsorbable solutes, the fraction of the water load passively reabsorbed in the proximal tubule is equal to the fraction of the sodium load reabsorbed (Uttamsingh et al., 1985). Taking into account details of laboratory measurements, this fraction also includes sodium reabsorbed in the loop of Henle (Seidlerová et al., 2006). The rate of fluid reabsorption from the distal tubules and collecting ducts depends on the influence of plasma vasopressin concentration, which was determined in the Uttamsingh model (Uttamsingh et al., 1985) on the basis of experimental measurements (Dehaven and Shapiro, 1970).

• Blood viscosity is calculated according to the study (Hund et al., 2017) as a function of hematocrit. We took into account that values of blood flows through the systemic and pulmonary microvessels should be inversely proportional to the normalized blood viscosity (Guyton et al., 1972).

Aldosterone Secretion

When modeling the effects of potassium and sodium on the secretion of aldosterone, we considered the following facts:

• Extremely low potassium levels (2 mmol/l) actually reduce aldosterone production stimulated by angiotensin II, but does not stop it completely (Kojima et al., 1985; Chen et al., 1999).

• A low-sodium diet does not directly affect the aldosterone secretion, but indirectly through activation of the RAAS, upregulation of AT1 receptor levels, and hyperplasia of the zona glomerulosa (Bollag, 2014).

• A rise of 1 mmol/l in serum potassium concentration doubles the aldosterone secretion (Pralong et al., 1992; Bollag, 2014).

Glomerular Dynamics

The total vascular resistance through the kidneys (RVR) is determined by the sum of the resistances in the individual vasculature segments (Hall, 2011). Table 1 shows the different variants of the formula for the RVR calculation used in the different studies. Note, that Gómez (1951) gives out the major renal resistances which can be estimated from clinically available data. Therefore, following this research, we used the same formula for RVR and applied it to derive an equation for renal blood flow (Supplementary Appendix B). We also included to the model such clinical variables as effective renal plasma flow, filtration fraction, and total protein (Škrtić et al., 2015).

TABLE 1

Table 1. The RVR determination as the sum of the resistances of separate vessels: A, interlobar, arcuate, and interlobular arteries; B, afferent arterioles; C, glomerular capillaries; D, efferent arterioles; E, peritubular capillaries; F, interlobar, interlobular, and arcuate veins.

Starting Values of Variables

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):

$${\mathit{\text{a}}}_{\mathit{i}}={\mathit{\text{c}}}_{\mathit{i}}\cdot m^{b_i}$$

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:

$$\begin{array}{c}{\mathit{\text{V}}}_{HL}=\mathbf{0.03}\cdot \mathit{\text{V}},{\mathit{\text{V}}}_{\mathit{AL}}=\mathbf{0.13}\cdot \mathit{\text{V}},\hfill \\ {\mathit{\text{V}}}_{\mathit{HR}}=\mathbf{0.03}\cdot \mathit{\text{V}},{\mathit{\text{V}}}_{\mathit{AR}}=\mathbf{0.035}\cdot \mathit{\text{V}}\hfill \end{array}$$

$$\begin{array}{c}{\mathit{\text{V}}}_{\mathit{VR}}=\mathbf{0.065}\cdot \mathit{\text{V}},\hfill \\ {\mathit{\text{V}}}_{\mathit{VL}}=\mathit{\text{V}}-V_{\mathit{AL}}-{\mathit{\text{V}}}_{\mathit{AR}}-V_{\mathit{HL}}-{\mathit{\text{V}}}_{\mathit{HR}}-{\mathit{\text{V}}}_{\mathit{VR}}\hfill \end{array}$$

The corresponding unstressed volumes we determined as:

$${\mathrm{\omega }}_{\mathit{i}}={\mathit{\text{k}}}_{\mathit{i}}\cdot {\mathit{\text{V}}}_{\mathit{i}},\mathit{\text{i}}\in \{\mathit{HL},\mathit{AL},\mathit{VL},\mathit{HR},\mathit{AR},\mathit{VR}\}$$

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).

The Model Calibration

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:

• Patient 1: Uncomplicated hypertension (Ferlinz, 1980).

• Patient 2: non-hypertensive diastolic heart failure (Fujimoto et al., 2008).

• Patient 3: Hypertension and LV hypertrophy without heart failure (Melenovsky et al., 2007).

• Patient 4: Hypertensive diastolic heart failure without LV hypertrophy (Fujimoto et al., 2008).

• Patient 5: Pulmonary hypertension and left heart disease (Wright et al., 2017).

The model calibration included two main steps.

(1) Search for the equilibrium state matching physiology values in a healthy human: We considered the normal value ranges for clinically measurable quantities used in the model. Such quantities cover 49 of 132 model parameters (Supplementary Appendix Table A3) and 69 of 160 model variables (Supplementary Appendix Table A4). The remaining parameters and variables of the model either cannot be measured in the laboratory, or we could not find available data to estimate their ranges. Then, 132 parameters of the model were determined in the following way.

(a) The values of 92 parameters were taken from the basic models (Guyton et al., 1972; Karaaslan et al., 2005; Proshin and Solodyannikov, 2006; Paeme et al., 2011; Hallow et al., 2014; Hallow and Gebremichael, 2017a; Rosalina et al., 2019) or found in the research articles (Gómez, 1951; Škrtić et al., 2015; Digne-Malcolm et al., 2016; Hund et al., 2017; Neal et al., 2018).

(b) The rest 7 new and 33 reused parameters, whose values directly affected the compliance with the normal ranges, were fitted by zeroing of the penalty function (5) using the constrained optimization evolutionary algorithm (see section “Materials and Methods”).

(2) Calibration of the model to the pathological equilibrium states.

(a) We defined a set of 65 fitting parameters, whose values can reasonably vary depending on the diseases of the abstract patients given above (Supplementary Appendix Table C1). These parameters are either directly related to the diseases, and therefore their values can deviate from the norm, or vary within the known normal ranges, if the diseases do not affect them.

(b) We determined a set of 57 physiological constraints imposed on the model variables in conformity with the clinical measurements in the diseases or in the normal state, if there was no available data confirming the effect of the diseases on these variables (Supplementary Appendix Table C2). Note that at the first step, we considered constraints for 69 clinical variables, whereas at this step the number of constraints is lower. This is due to the fact that we eliminated some of the excessive constraints to accelerate the optimization process. For example, the heart rate is related to the cardiac cycle length. Thus, we took into account the physiological range for the first variable and excluded the range for the second one. The cardiac output is equal to the product of the stroke volume and the heart rate. Therefore, we considered the constraints for the first two variables and skipped the constraints for the last one, etc.

(c) We chose standard patient characteristics provided in clinical studies of the diseases (references in definition of patients 1–5): age, weight, height, gender, systolic/diastolic blood pressure, heart rate, stroke volume, ejection fraction, and hemoglobin. We used the average values of these characteristics reported in the studies (upon availability): weight and hemoglobin were directly set to the model; weight, height, and gender were taken to evaluate total body water of abstract patients (by the Nadler equation, Nadler et al., 1962); remaining values were used to estimate the corresponding model variables by the minimization of the distance function (4).

(d) We specified the optimization problem with the set of fitting parameters (a), the penalty function (b), and the objective function (c), and solved this problem as described in the section “Materials and Methods” for each test case.

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).

TABLE 2

Table 2. Equilibrium values representing abstract patients with different cardiovascular diseases.

Use Cases

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).

Comparison of Different Physiological States

As the typical examples, consider the following:

• Comparison of normal and abnormal physiology (uncomplicated hypertension vs. normal state).

• Comparison of states with the same disease but different underlying pathophysiologic mechanisms (non-hypertensive vs. hypertensive diastolic heart failure).

• Comparison of the pathological states reproducing complex cardiovascular diseases (pulmonary hypertension and left heart disease vs. systemic hypertension and left ventricular hypertrophy).

Figure 6 shows the sample plots of the left ventricular pressure-volume loops simulated for all these cases. Similar plots can be automatically generated by the code in JavaScript included into the Jupyter file for any equilibrium states of the model and any variables of interest. Specifically, in the given plots, we can see the dynamics corresponding to the considered patient diseases:

• In the untreated hypertensive patients, left ventricular systolic pressure can be substantially higher than in normal subjects (Antony et al., 1993).

• Left ventricular end-systolic pressure is higher in hypertensive than in non-hypertensive patients with diastolic heart failure, while left ventricular end-diastolic pressure is elevated in the both groups (Fujimoto et al., 2008).

• Patients with congestive heart failure (generally accompanying pulmonary hypertension) can have significantly higher left ventricular end-diastolic and end-systolic volumes (Mehta et al., 2000).

FIGURE 6

Figure 6. Comparison of left ventricular pressure-volume loops simulated for different equilibrium states of the model. (A) Uncomplicated hypertension vs. normal state. (B) Non-hypertensive vs. hypertensive diastolic heart failure (DHF). (C) Pulmonary hypertension and left heart disease (PH/LHD) vs. systemic hypertension and left ventricular hypertrophy (SH/LVH).

Comparison of Sensitivity of Different Patients to the Change in Physiological Conditions

As the test case of this problem, we used the experiment with varying sodium intake presented in He et al. (2001). The study by He et al. (2001) involved a 5-day high sodium diet (≈ 350 mmol/d) followed by a 5-day low sodium diet (10–20 mmol/d) in normotensive and hypertensive individuals.

We designed the corresponding simulation experiment by the following way:

• 0–1 day – the normal diet with initial equilibrium values of sodium intake in the model states.

• 1–6 days – high sodium diet with sodium intake equal to 0.24306 mEq/min ≈ 350 mmol/d.

• 6–11 days – low sodium diet with sodium intake equal to 0.01042 mEq/min ≈ 15 mmol/d.

• 11–21 days – the normal diet returning the model dynamics to the initial equilibrium.

Figure 7 shows the dynamic results of the model in all states for relative values (normalized to the initial value) of sodium intake, plasma sodium, mean arterial pressure, heart rate, plasma renin activity and aldosterone concentration. Table 3 reveals the exact and relative difference in these values between high (6 days) and low (11 days) salt diets. As can be seen from these data, hypertensive patients have a greater fall in blood pressure from high-salt to low-salt diet than normotensive subjects. This result is consistent with the conclusions by He et al. (2001). However, the laboratory measurements by the authors also demonstrated that supine pulse rate on average did not change with acute salt restriction, whereas plasma renin activity and aldosterone concentration had higher growth in normotensive than in hypertensive individuals. In our simulation test, we observe another dynamic: a decrease in blood pressure is accompanied by an increase in heart rate. Such dynamics is in line with the fact that sodium reduction can increase heart rate (Graudal et al., 2016). For plasma renin activity and aldosterone concentration, we get the results opposite to conclusions by He et al. (2001). The model shows higher growth of these quantities in hypertensive patients. However, while plasma renin activity responses on average are stronger in normotensive than in hypertensive populations (Graudal et al., 2017), it can be weaker in some individual cases (Graudal et al., 2021). Therefore, the generated model states give acceptable dynamics.

FIGURE 7

Figure 7. Simulation results of the sodium load experiment for all states of the model. The simulation experiment was designed in accordance with experimental study by He et al. (2001) and involves a 5-day high sodium diet followed by a 5-day low sodium diet. Variable values are normalized to the initial equilibrium values declared in the states.

TABLE 3

Table 3. Exact and relative difference in sodium intake (mEq/min), plasma sodium (mEq/l), mean arterial pressure (MAP, mmHg), heart rate (HR, beats/min), plasma renin activity (PRA, fmol/ml/min) and plasma aldosterone concentration (PAC, pg/ml) between high (day 6) and low (day 11) salt diets.

Discussion

The main goal of our investigation is the modeling of a hypertensive human personal response to antihypertensive therapy. Since arterial hypertension is the multifactorial disease, which can take different forms (pulmonary or systemic) and can be complicated by related pathologies, such as heart or renal failure, the modeling process is divided into several steps.

• Step 1: Creation of a mathematical model of human physiology with a level of detail sufficient to study the issues of arterial hypertension. True personalization of drug therapies should rely on a virtual patient, the digital twin of a real individual, which is formed and accumulated throughout his life as a result of interaction with the health care system (Lehrach, 2016). We suppose that it is unrealistic to build a virtual patient for all occasions now. Thus, our approach is to construct a set of basic modules (blocks) and assemble a model from them (as from Lego blocks) for a given patient and disease.

• Step 2: Pharmacokinetic/pharmacodynamic (PK/PD) modeling of antihypertensive drugs. First-line antihypertensive medications include angiotensin-converting enzyme inhibitors, angiotensin II receptor blockers, calcium-channel blockers, thiazide diuretics and β-Adrenoreceptor blockers (Oparil et al., 2018). It is possible to determine the certain points of influence on the model for each of these classes, construct the corresponding PK/PD models and validate unknown dynamic constants using related clinical trials (Hallow et al., 2014).

• Step 3: Personalizing the model. Each equilibrium parameterization of the model within physiological ranges can be considered as a virtual patient. To relate him with a real person, we can use some values from the laboratory analyses. However, this allows us to get only a small part of the model quantities. To solve the problem with unknown personal parameters, we can build a set of virtual patients and consider significant variation of unknown physiological values. Treatment simulation of such a virtual population makes it possible to identify virtual groups with a similar reaction to the drugs and analyze which features of the patient can contribute to the effectiveness (or ineffectiveness) of the antihypertensive therapy.

This article provides implementation of the first step and presents the model of cardiovascular and renal systems. As the basis for the model construction, we used early published models developed for each of these physiology systems separately (Karaaslan et al., 2005, 2014; Proshin and Solodyannikov, 2006; Hallow et al., 2014; Hallow and Gebremichael, 2017a). Note that the reduction in dietary salt intake leads to a decrease in blood pressure and therefore, can be viewed as a simple antihypertensive therapy (Mahtani, 2009; Frisoli et al., 2012). As follows from Figure 7, the model reproduces the dynamics consistent with this statement. As for the modeling of the individual response to complex therapy including antihypertensive medications with different mechanisms of action (Step 2 and 3), this is a task for future studies.

The Model Strength

Two Physiological Systems Instead of One

The strong connection between renal and cardiovascular disease reflects the complex interactions between heart and kidneys (Stefanadis, 2010). Therefore, it is reasonable to model in detail these body systems together. In relation to the study of personal response to antihypertensive therapies, simultaneous consideration of these systems allows simulating virtual patients with different combinations of heart and renal diseases associated with hypertension. This also provides a tool for modeling the hypotensive effect in response to β-blockers treatment, which is not possible with consideration of the renal function alone (Hallow et al., 2014).

Modular Representation

The modular approach facilitates development of complex models by representing them as combinations of submodels (Kiselev and Kolpakov, 2013). The structure of the model is clear and understandable. The model can be easily expanded. The separate modules can be independently modified and improved.

Agent-Based Approach

Agent-based models are unique in their ability to integrate combinations of heterogeneous processes and investigate their respective dynamics. These models are flexible in their execution and permit the aggregation of processes across time scales (Glen et al., 2019). So in our case, the cardiovascular system modules are simulated in fraction of seconds which is due to the work of the heart, while the renal system modules are measured in minutes.

Wide Range of Parameter Values

Our model is not geared toward parametrizations representing a limited group of healthy people or patients with the same combination of diseases, but has a wide parameter space, which is sufficient for simulation of patients with systemic/pulmonary arterial hypertension and main concomitant cardiovascular diseases.

The Model Limitations

Limitations by the Parameter Constraints

Parameter bounds and constraints, which we collected on the basis of known clinical studies and listed in the Supplementary Appendix Tables C1, C2, are suitable for simulation of a wide range of patients with cardiovascular diseases, but still do not reflect all possible cases and require further extension and systematization to detail physiological ranges depending on patient features.

Limitations by the Model Equations

We use the same form of equations to simulate different pathological conditions. However, in some diseases, a part of variables may be outside the acceptable ranges. For example, in the study (Riegger et al., 1982), one group of patients with congestive heart failure had inappropriately high values of plasma antidiuretic hormone (14.5 ± 8.8 pg/ml) in relation to their plasma osmolality, which was well below normal values (276 ± 23 mOsmol/kg water). At the same time, the equation used in the model, does not allow the analysis of osmolality values less than 271 mOsmol/kg. Thus, to simulate such patients, the function of calculation of plasma antidiuretic hormone must be advanced.

Possible Directions of the Model Development

Aging as a Key Factor in Cardiovascular Diseases

Aging-mediated structural and biochemical modifications coupled with gradual loss of autonomic nervous system regulation and vascular stiffening are consistently implicated in the progressive increase in mechanical burden and functional breakdown of the heart and vessels (Fajemiroye et al., 2018). Thus, it seems to us very important to introduce age in the equations modeling the age-dependent variables.

Genetic Contribution

Estimated heritability of systolic and diastolic blood pressure lies in the ranges of 15–40% and 15–30%, respectively (Delles and Padmanabhan, 2012). More than 900 known genetic loci indicate that sites for blood pressure control involve various organs, including the kidneys and nervous system (Lin et al., 2020). Therefore, extension of the model to take into account personal data of genetic testing can help resolve the questions regarding individual blood pressure regulation.

Personalized Cardiac Electrophysiology Modeling

To personalize the cardiac parameters of the model, it is possible to use data obtained from clinical imaging. However, when imaging data is noisy, an alternative rule-based methodology can be utilized to simulate electrical wave propagation and mechanical contraction in the heart (Bayer et al., 2012; Lopez-Perez et al., 2015; Nguyen et al., 2020). Application of this methodology can be useful for extension the model to involve the personalized anatomy of the heart.

Finally, note that we did not aim to reproduce all possible variations of cardiovascular diseases. Nevertheless, we suppose that the model can be easily extended to any group of patients, depending on the study purpose. Summarizing the above, we believe that the composite model of cardiovascular and renal systems could be useful for further investigation of cardiovascular diseases and drug development. The BioUML implementation of the model is available at: https://gitlab.sirius-web.org/virtual-patient/blood-pressure-regulation.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: gitlab.sirius-web.org/virtual-patient/blood-pressure-regulation.

Author Contributions

IK and EK designed the agent-based model of cardiovascular and renal systems, and implemented necessary tools in the BioUML software (agent-based modeling and simulation tools, and optimization tools). EK validated the model based on the clinical data from scientific literature. RS consulted on the physiological questions. GL consulted on the questions of hypertension and antihypertensive therapy, coordinated the medical part of the study. FK coordinated creation of the model and development of BioUML. EK and RS prepared of the model documentation (Supplementary Material). EK, IK, RS, GL, and FK wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was funded in part by the Sirius University of Science and Technology.

Conflict of Interest

EK, IK, RS, and FK were employed by company Biosoft.Ru, Ltd., which develops and supports the BioUML platform. During the research the BioUML platform was improved by authors for the purpose of the research.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors are deeply grateful to V. I. Baranov, who made a great contribution to the model development and discussion and became a victim of COVID-19.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2021.746300/full#supplementary-material

Footnotes

1. ^ www.samara-dialog.ru/help/eng/help.htm
2. ^ www.biouml.org
3. ^ jupyter.org

References

Abram, S., Hodnett, B. L., Summers, R. L., Coleman, T. G., and Hester, R. L. (2007). Quantitative circulatory physiology: an integrative mathematical model of human physiology for medical education. Adv. Physiol. Educ. 31, 202–210. doi: 10.1152/advan.00114.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, A. M., Zhuo, J., and Mendelsohn, F. A. (2000). Localization and function of angiotensin AT1 receptors. Am. J. Hypertens. 13, 31S–38S. doi: 10.1016/s0895-7061(99)00249-6

CrossRef Full Text | Google Scholar

Alon, U. (2003). Biological networks: the tinkerer as an engineer. Science 301, 1866–1867. doi: 10.1126/science.1089072

PubMed Abstract | CrossRef Full Text | Google Scholar

Antony, I., Nitenberg, A., Foult, J. M., and Aptecar, E. (1993). Coronary vasodilator reserve in untreated and treated hypertensive patients with and without left ventricular hypertrophy. J. Am. Coll. Cardiol. 22, 514–520. doi: 10.1016/0735-1097(93)90058-9

CrossRef Full Text | Google Scholar

Ardaillou, R., Chansel, D., Chatziantoniou, C., and Dussaule, J. C. (1999). Mesangial AT1 receptors: expression, signaling, and regulation. J. Am. Soc. Nephrol. Suppl. 11, S40–S46.

Google Scholar

Bandulik, S., Tauber, P., Lalli, E., Barhanin, J., and Warth, R. (2015). Two-pore domain potassium channels in the adrenal cortex. Pflugers Arch. Eur. J. Physiol. 467, 1027–1042. doi: 10.1007/s00424-014-1628-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bayer, J. D., Blake, R. C., Plank, G., and Trayanova, N. A. (2012). A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. Ann. Biomed. Eng. 40, 2243–2254. doi: 10.1007/s10439-012-0593-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Beard, D. A. (2018). Assessing the validity and utility of the Guyton model of arterial blood pressure control. Hypertension 72, 1272–1273. doi: 10.1161/HYPERTENSIONAHA.118.11998

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhagat, C. I., Garcia-Webb, P., Fletcher, E., and Beilby, J. P. (1984). Calculated vs measured plasma osmolalities revisited. Clin. Chem. 30, 1703–1705.

Google Scholar

Blantz, R. C., Konnen, K. S., and Tucker, B. J. (1976). Angiotensin II effects upon the glomerular microcirculation and ultrafiltration coefficient of the rat. J. Clin. Invest. 57, 419–434. doi: 10.1172/JCI108293

PubMed Abstract | CrossRef Full Text | Google Scholar

Blinov, M. L., Ruebenacker, O., and Moraru, I. I. (2008). Complexity and modularity of intracellular networks: a systematic approach for modelling and simulation. IET Syst. Biol. 2, 363–368. doi: 10.1049/iet-syb:20080092

CrossRef Full Text | Google Scholar

Bollag, W. B. (2014). Regulation of aldosterone synthesis and secretion. Compr. Physiol. 4, 1017–1055. doi: 10.1002/cphy.c130037

PubMed Abstract | CrossRef Full Text | Google Scholar

Boogers, M. J., van Werkhoven, J. M., Schuijf, J. D., Delgado, V., El-Naggar, H. M., Boersma, E., et al. (2011). Feasibility of diastolic function assessment with cardiac CT: feasibility study in comparison with tissue Doppler imaging. JACC Cardiovasc. Imaging 4, 246–256. doi: 10.1016/j.jcmg.2010.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Brasch, H., Sieroslawski, L., and Dominiak, P. (1993). Angiotensin II increases norepinephrine release from atria by acting on angiotensin subtype 1 receptors. Hypertension 22, 699–704. doi: 10.1161/01.hyp.22.5.699

CrossRef Full Text | Google Scholar

Brown, P. N., Byrne, G. D., and Hindmarsh, A. C. (1989). VODE: a variable-coefficient ODE solver. SIAM J. Sci. Stat Comput. 10, 1038–1051. doi: 10.1137/0910062

CrossRef Full Text | Google Scholar

Caudron, J., Fares, J., Bauer, F., and Dacher, J. N. (2011). Evaluation of left ventricular diastolic function with cardiac MR imaging. Radiographics 31, 239–261. doi: 10.1148/rg.311105049

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X. L., Bayliss, D. A., Fern, R. J., and Barrett, P. Q. (1999). A role for T-type Ca²⁺ channels in the synergistic control of aldosterone production by ANG II and K⁺. Am. J. Physiol. 276, F674–F683. doi: 10.1152/ajprenal.1999.276.5.F674

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, Y., Thalhauser, C. J., Smithline, S., Pagidala, J., Miladinov, M., Vezina, H. E., et al. (2017). QSP toolbox: computational implementation of integrated workflow components for deploying multi-scale mechanistic models. AAPS J. 19, 1002–1016. doi: 10.1208/s12248-017-0100-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, S. D., and Hindmarsh, A. C. (1996). CVODE, a stiff/nonstiff ODE solver in C. Comput. Phys. 10, 138–143. doi: 10.1063/1.4822377

CrossRef Full Text | Google Scholar

Dehaven, J. C., and Shapiro, N. Z. (1970). Simulation of the renal effects of antidiuretic hormone (ADH) in man. J. Theor. Biol. 28, 261–286. doi: 10.1016/0022-5193(70)90055-x

CrossRef Full Text | Google Scholar

Delles, C., and Padmanabhan, S. P. (2012). Genetics and hypertension: is it time to change my practice? Can. J. Cardiol. 28, 296–304. doi: 10.1016/j.cjca.2012.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Digne-Malcolm, H., Frise, M. C., and Dorrington, K. L. (2016). How do antihypertensive drugs work? Insights from studies of the renal regulation of arterial blood pressure. Front. Physiol. 7:320. doi: 10.3389/fphys.2016.00320

PubMed Abstract | CrossRef Full Text | Google Scholar

Dorwart, W. V., and Chalmers, L. (1975). Comparison of methods for calculating serum osmolality form chemical concentrations, and the prognostic value of such calculations. Clin. Chem. 21, 190–194.

Google Scholar

Edelman, I. S., Leibman, J., O’meara, M. P., and Birkenfeld, L. W. (1958). Interrelations between serum sodium concentration, serum osmolarity and total exchangeable sodium, total exchangeable potassium and total body water. J. Clin. Invest. 37, 1236–1256. doi: 10.1172/JCI103712

PubMed Abstract | CrossRef Full Text | Google Scholar

Fajemiroye, J. O., da Cunha, L. C., Saavedra-Rodríguez, R., Rodrigues, K. L., Naves, L. M., Mourão, A. A., et al. (2018). Aging-induced biological changes and cardiovascular diseases. Biomed. Res. Int. 2018:7156435. doi: 10.1155/2018/7156435

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferlinz, J. (1980). Right ventricular performance in essential hypertension. Circulation 61, 156–162. doi: 10.1161/01.cir.61.1.156

CrossRef Full Text | Google Scholar

Frisoli, T. M., Schmieder, R. E., Grodzicki, T., and Messerli, F. H. (2012). Salt and hypertension: is salt dietary reduction worth the effort? Am. J. Med. 125, 433–439. doi: 10.1016/j.amjmed.2011.10.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujimoto, N., Onishi, K., Dohi, K., Tanabe, M., Kurita, T., Takamura, T., et al. (2008). Hemodynamic characteristics of patients with diastolic heart failure and hypertension. Hypertens. Res. 31, 1727–1735. doi: 10.1291/hypres.31.1727

PubMed Abstract | CrossRef Full Text | Google Scholar

Gan, C., Lankhaar, J. W., Marcus, J. T., Westerhof, N., Marques, K. M., Bronzwaer, J. G., et al. (2006). Impaired left ventricular filling due to right-to-left ventricular interaction in patients with pulmonary arterial hypertension. Am. J. Physiol. Heart Circ. Physiol. 290, H1528–H1533. doi: 10.1152/ajpheart.01031.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Gazioglu, K., and Yu, P. N. (1967). Pulmonary blood volume and pulmonary capillary blood volume in valvular heart disease. Circulation 35, 701–709. doi: 10.1161/01.cir.35.4.701

CrossRef Full Text | Google Scholar

Glen, C. M., Kemp, M. L., and Voit, E. O. (2019). Agent-based modeling of morphogenetic systems: advantages and challenges. PLoS Comput. Biol. 15:e1006577. doi: 10.1371/journal.pcbi.1006577

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez, D. M. (1951). Evaluation of renal resistances, with special reference to changes in essential hypertension. J. Clin. Invest. 30, 1143–1155. doi: 10.1172/JCI102534

PubMed Abstract | CrossRef Full Text | Google Scholar

González-Abraldes, J., Albillos, A., Bañares, R., Del Arbol, L. R., Moitinho, E., Rodríguez, C., et al. (2001). Randomized comparison of long-term losartan versus propranolol in lowering portal pressure in cirrhosis. Gastroenterology 121, 382–388. doi: 10.1053/gast.2001.26288

PubMed Abstract | CrossRef Full Text | Google Scholar

Graudal, N., Hubeck-Graudal, T., and Jurgens, G. (2021). Influence of sodium intake and change in sodium intake on plasma-renin in man. EClinicalMedicine 33:100750. doi: 10.1016/j.eclinm.2021.100750

PubMed Abstract | CrossRef Full Text | Google Scholar

Graudal, N. A., Hubeck-Graudal, T., and Jürgens, G. (2016). Reduced dietary sodium intake increases heart rate. a meta-analysis of 63 randomized controlled trials including 72 study populations. Front. Physiol. 7:111. doi: 10.3389/fphys.2016.00111

PubMed Abstract | CrossRef Full Text | Google Scholar

Graudal, N. A., Hubeck-Graudal, T., and Jurgens, G. (2017). Effects of low sodium diet versus high sodium diet on blood pressure, renin, aldosterone, catecholamines, cholesterol, and triglyceride. Cochrane Database Syst. Rev. 4:CD004022. doi: 10.1002/14651858.CD004022.pub3

PubMed Abstract | CrossRef Full Text | Google Scholar

Guyton, A. C., Coleman, T. G., and Granger, H. J. (1972). Circulation: overall regulation. Annu. Rev. Physiol. 34, 13–46. doi: 10.1146/annurev.ph.34.030172.000305

PubMed Abstract | CrossRef Full Text | Google Scholar

Hall, J. E. (2011). Guyton and Hall textbook of Medical Physiology, 12th Edn. Philadelphia, PA: Saunders Elsevier, 1091.

Google Scholar

Hallow, K. M., and Gebremichael, Y. (2017a). A quantitative systems physiology model of renal function and blood pressure regulation: model description. CPT Pharmacometrics Syst. Pharmacol. 6, 383–392. doi: 10.1002/psp4.12178

PubMed Abstract | CrossRef Full Text | Google Scholar

Hallow, K. M., and Gebremichael, Y. (2017b). A quantitative systems physiology model of renal function and blood pressure regulation: application in salt-sensitive hypertension. CPT Pharmacometrics Syst. Pharmacol. 6, 393–400. doi: 10.1002/psp4.12177

PubMed Abstract | CrossRef Full Text | Google Scholar

Hallow, K. M., Lo, A., Beh, J., Rodrigo, M., Ermakov, S., Friedman, S., et al. (2014). A model-based approach to investigating the pathophysiological mechanisms of hypertension and response to antihypertensive therapies: extending the Guyton model. Am. J. Physiol. Regul. Integr. Comp. Physiol. 306, R647–R662. doi: 10.1152/ajpregu.00039.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Hallow, K. M., Van Brackle, C. H., Anjum, S., and Ermakov, S. (2021). Cardiorenal systems modeling: left ventricular hypertrophy and differential effects of antihypertensive therapies on hypertrophy regression. Front. Physiol. 12:679930. doi: 10.3389/fphys.2021.679930

PubMed Abstract | CrossRef Full Text | Google Scholar

Hammer, M., Ladefoged, J., and Olgaard, K. (1980). Relationship between plasma osmolality and plasma vasopressin in human subjects. Am. J. Physiol. 238, E313–E317. doi: 10.1152/ajpendo.1980.238.4.E313

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartwell, L. H., Hopfield, J. J., Leibler, S., and Murray, A. W. (1999). From molecular to modular cell biology. Nature 402(Suppl. 6761), C47–C52. doi: 10.1038/35011540

PubMed Abstract | CrossRef Full Text | Google Scholar

He, F. J., Markandu, N. D., and MacGregor, G. A. (2001). Importance of the renin system for determining blood pressure fall with acute salt restriction in hypertensive and normotensive whites. Hypertension. 38, 321–325. doi: 10.1161/01.hyp.38.3.321

CrossRef Full Text | Google Scholar

Hernández, A. I., Le Rolle, V., Defontaine, A., and Carrault, G. (2009). A multiformalism and multiresolution modelling environment: application to the cardiovascular system and its regulation. Philos. Trans. A Math. Phys. Eng. Sci. 367, 4923–4940. doi: 10.1098/rsta.2009.0163

PubMed Abstract | CrossRef Full Text | Google Scholar

Hester, R. L., Brown, A. J., Husband, L., Iliescu, R., Pruett, D., Summers, R., et al. (2011). HumMod: a modeling environment for the simulation of integrative human physiology. Front. Physiol. 2:12. doi: 10.3389/fphys.2011.00012

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirschfeld, S., Meyer, R., Schwartz, D. C., Korfhagen, J., and Kaplan, S. (1975). Measurement of right and left ventricular systolic time intervals by echocardiography. Circulation 51, 304–309. doi: 10.1161/01.cir.51.2.304

CrossRef Full Text | Google Scholar

Hoops, S., Sahle, S., Gauges, R., Lee, C., Pahle, J., Simus, N., et al. (2006). COPASI — a complex pathway simulator. Bioinformatics 22, 3067–3074. doi: 10.1093/bioinformatics/btl485

PubMed Abstract | CrossRef Full Text | Google Scholar

Hucka, M., Bergmann, F. T., Chaouiya, C., Dräger, A., Hoops, S., Keating, S. M., et al. (2019). The systems biology markup language (SBML): language specification for level 3 version 2 core. J. Integr. Bioinform. 16:20190021. doi: 10.1515/jib-2017-0081

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughes, A. D. (1998). Molecular and cellular mechanisms of action of angiotensin II (AT1) receptors in vascular smooth muscle. J. Hum. Hypertens. 12, 275–281. doi: 10.1038/sj.jhh.1000635

PubMed Abstract | CrossRef Full Text | Google Scholar

Hund, S. J., Kameneva, M. V., and Antaki, J. F. (2017). A quasi-mechanistic mathematical representation for blood viscosity. Fluids 2:10. doi: 10.3390/fluids2010010

CrossRef Full Text | Google Scholar

Ikeda, N., Marumo, F., Shirataka, M., and Sato, T. (1979). A model of overall regulation of body fluids. Ann. Biomed. Eng. 7, 135–166. doi: 10.1007/BF02363132

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamishirado, H., Inoue, T., Fujito, T., Kase, M., Shimizu, M., Sakai, Y., et al. (1997). Effect of enalapril maleate on cerebral blood flow in patients with chronic heart failure. Angiology 48, 707–713. doi: 10.1177/000331979704800806

PubMed Abstract | CrossRef Full Text | Google Scholar

Karaaslan, F., Denizhan, Y., and Hester, R. (2014). A mathematical model of long-term renal sympathetic nerve activity inhibition during an increase in sodium intake. Am. J. Physiol. Regul. Integr. Comp. Physiol. 306, R234–R247. doi: 10.1152/ajpregu.00302.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Karaaslan, F., Denizhan, Y., Kayserilioglu, A., and Gulcur, H. O. (2005). Long-term mathematical model involving renal sympathetic nerve activity, arterial pressure, and sodium excretion. Ann. Biomed. Eng. 33, 1607–1630. doi: 10.1007/s10439-005-5976-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kassab, G., and Guccione, J. (2019). Editorial: mathematical modeling of cardiovascular systems: from physiology to the clinic. Front. Physiol. 10:1259. doi: 10.3389/fphys.2019.01259

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiselev, I., and Kolpakov, F. (2013). Modular modeling of biological systems. Virtual Biol. 1:e11. doi: 10.12704/vb/e11

CrossRef Full Text | Google Scholar

Klinke, D. J. II (2015). Enhancing the discovery and development of immunotherapies for cancer using quantitative and systems pharmacology: interleukin-12 as a case study. J. Immunother. Cancer 3:27. doi: 10.1186/s40425-015-0069-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, M., Furukawa, Y., and Chiba, S. (1978). Positive chronotropic and inotropic effects of angiotensin II in the dog heart. Eur. J. Pharmacol. 50, 17–25. doi: 10.1016/0014-2999(78)90249-2

CrossRef Full Text | Google Scholar

Kojima, I., Kojima, K., and Rasmussen, H. (1985). Role of calcium fluxes in the sustained phase of angiotensin II-mediated aldosterone secretion from adrenal glomerulosa cells. J. Biol. Chem. 260, 9177–9184.

Google Scholar

Kolpakov, F., Akberdin, I., Kashapov, T., Kiselev, I., Kolmykov, S., Kondrakhin, Y., et al. (2019). BioUML: an integrated environment for systems biology and collaborative analysis of biomedical data. Nucleic Acids Res. 47, W225–W233. doi: 10.1093/nar/gkz440

PubMed Abstract | CrossRef Full Text | Google Scholar

Konstam, M. A., Kronenberg, M. W., Rousseau, M. F., Udelson, J. E., Melin, J., Stewart, D., et al. (1993). Effects of the angiotensin converting enzyme inhibitor enalapril on the long-term progression of left ventricular dilatation in patients with asymptomatic systolic dysfunction. SOLVD (Studies of Left Ventricular Dysfunction) Investigators. Circulation 88(part 1), 2277–2283. doi: 10.1161/01.cir.88.5.2277

CrossRef Full Text | Google Scholar

Kutumova, E. O., Kiselev, I. N., Sharipov, R. N., Lavrik, I. N., and Kolpakov, F. A. (2012). A modular model of the apoptosis machinery. Adv. Exp. Med. Biol. 736, 235–245. doi: 10.1007/978-1-4419-7210-1_13

CrossRef Full Text | Google Scholar

Le Novère, N., Hucka, M., Mi, H., Moodie, S., Schreiber, F., Sorokin, A., et al. (2009). The systems biology graphical notation. Nat. Biotechnol. 27, 735–741. doi: 10.1038/nbt.1558

PubMed Abstract | CrossRef Full Text | Google Scholar

Lehrach, H. (2016). Omics approaches to individual variation: modeling networks and the virtual patient. Dialogues Clin. Neurosci. 18, 253–265. doi: 10.31887/DCNS.2016.18.3/hlehrach

CrossRef Full Text | Google Scholar

Leighton, R. F., Weissler, A. M., Weinstein, P. B., and Wooley, C. F. (1971). Right and left ventricular systolic time intervals. Effects of heart rate, respiration and atrial pacing. Am. J. Cardiol. 27, 66–72. doi: 10.1016/0002-9149(71)90084-1

CrossRef Full Text | Google Scholar

Lin, H. J., Guo, X., and Rotter, J. I. (2020). “6 - The genetics of blood pressure regulation,” in Emery and Rimoin’s Principles and Practice of Medical Genetics and Genomics: Cardiovascular, Respiratory, and Gastrointestinal Disorders, 7th Edn, eds R. E. Pyeritz, B. R. Korf, and W. W. Grody (Amsterdam: Elsevier Inc), 577.

Google Scholar

Lopez-Perez, A., Sebastian, R., and Ferrero, J. M. (2015). Three-dimensional cardiac computational modelling: methods, features and applications. Biomed. Eng. Online. 14:35. doi: 10.1186/s12938-015-0033-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Macal, C. M., and North, M. J. (2010). Tutorial on agent-based modelling and simulation. J. Simul. 4, 151–162. doi: 10.1057/jos.2010.3

CrossRef Full Text | Google Scholar

Maceira, A. M., Prasad, S. K., Khan, M., and Pennell, D. J. (2006a). Normalized left ventricular systolic and diastolic function by steady state free precession cardiovascular magnetic resonance. J. Cardiovasc. Magn. Reson. 8, 417–426. doi: 10.1080/10976640600572889

PubMed Abstract | CrossRef Full Text | Google Scholar

Maceira, A. M., Prasad, S. K., Khan, M., and Pennell, D. J. (2006b). Reference right ventricular systolic and diastolic function normalized to age, gender and body surface area from steady-state free precession cardiovascular magnetic resonance. Eur. Heart J. 27, 2879–2888. doi: 10.1093/eurheartj/ehl336

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahtani, K. R. (2009). Simple advice to reduce salt intake. Br. J. Gen. Pract. 59, 786–787. doi: 10.3399/bjgp09X472755

PubMed Abstract | CrossRef Full Text | Google Scholar

Maxwell, G., and Mackay, C. (2008). Application of a systems biology approach to skin allergy risk assessment. Altern. Lab. Anim. 36, 521–556. doi: 10.1177/026119290803600510

PubMed Abstract | CrossRef Full Text | Google Scholar

McDaniel, M., Carter, J., Keller, J. M., White, S. A., and Baird, A. (2019). Open source pharmacokinetic/pharmacodynamic framework: tutorial on the BioGears engine. CPT Pharmacometrics Syst. Pharmacol. 8, 12–25. doi: 10.1002/psp4.12371

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehta, S., Liu, P. P., Fitzgerald, F. S., Allidina, Y. K., and Douglas Bradley, T. (2000). Effects of continuous positive airway pressure on cardiac volumes in patients with ischemic and dilated cardiomyopathy. Am. J. Respir. Crit. Care Med. 161, 128–134. doi: 10.1164/ajrccm.161.1.9903055

PubMed Abstract | CrossRef Full Text | Google Scholar

Melenovsky, V., Borlaug, B. A., Rosen, B., Hay, I., Ferruci, L., Morell, C. H., et al. (2007). Cardiovascular features of heart failure with preserved ejection fraction versus nonfailing hypertensive left ventricular hypertrophy in the urban Baltimore community: the role of atrial remodeling/dysfunction. J. Am. Coll. Cardiol. 49, 198–207. doi: 10.1016/j.jacc.2006.08.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Montani, J. P., and Vliet, B. N. V. (2009). Understanding the contribution of Guyton’s large circulatory model to long-term control of arterial pressure. Exp. Physiol. 94, 382–388. doi: 10.1113/expphysiol.2008.043299

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, F. D. (1967). Body composition and its measurement in vivo. Br. J. Surg. 54, 431–435. doi: 10.1002/bjs.1800541309

PubMed Abstract | CrossRef Full Text | Google Scholar

Moss, R., Grosse, T., Marchant, I., Lassau, N., Gueyffier, F., and Thomas, S. R. (2012). Virtual patients and sensitivity analysis of the Guyton model of blood pressure regulation: towards individualized models of whole-body physiology. PLoS Comput. Biol. 8:e1002571. doi: 10.1371/journal.pcbi.1002571

PubMed Abstract | CrossRef Full Text | Google Scholar

Murch, S. D., La Gerche, A., Roberts, T. J., Prior, D. L., MacIsaac, A. I., and Burns, A. T. (2015). Abnormal right ventricular relaxation in pulmonary hypertension. Pulm. Circ. 5, 370–375. doi: 10.1086/681268

PubMed Abstract | CrossRef Full Text | Google Scholar

Nadler, S. B., Hidalgo, J. U., and Bloch, T. (1962). Prediction of blood volume in normal human adults. Surgery 51, 224–232.

Google Scholar

Natarajan, A. R., Eisner, G. M., Armando, I., Browning, S., Pezzullo, J. C., Rhee, L., et al. (2016). The renin-angiotensin and renal dopaminergic systems interact in normotensive humans. J. Am. Soc. Nephrol. 27, 265–279. doi: 10.1681/ASN.2014100958

PubMed Abstract | CrossRef Full Text | Google Scholar

Neal, C. R., Arkill, K. P., Bell, J. S., Betteridge, K. B., Bates, D. O., Winlove, C. P., et al. (2018). Novel hemodynamic structures in the human glomerulus. Am. J. Physiol. Renal Physiol. 315, F1370–F1384. doi: 10.1152/ajprenal.00566.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Neal, M. L., Cooling, M. T., Smith, L. P., Thompson, C. T., Sauro, H. M., Carlson, B. E., et al. (2014). A reappraisal of how to build modular, reusable models of biological systems. PLoS Comput. Biol. 10:e1003849. doi: 10.1371/journal.pcbi.1003849

PubMed Abstract | CrossRef Full Text | Google Scholar

Neubauer, B., Schrankl, J., Steppan, D., Neubauer, K., Sequeira-Lopez, M. L., Pan, L., et al. (2018). Ang II short loop feedback – is there a role of Ang II for the regulation of the renin system in vivo? Hypertension 71, 1075–1082. doi: 10.1161/HYPERTENSIONAHA.117.10357

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, M. K., and Kurtz, I. (2003). Are the total exchangeable sodium, total exchangeable potassium and total body water the only determinants of the plasma water sodium concentration? Nephrol Dial Transplant. 18, 1266–1271. doi: 10.1093/ndt/gfg112

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, T. D., Kadri, O. E., and Voronov, R. S. (2020). An introductory overview of image-based computational modeling in personalized cardiovascular medicine. Front. Bioeng. Biotechnol. 8:529365. doi: 10.3389/fbioe.2020.529365

PubMed Abstract | CrossRef Full Text | Google Scholar

Oparil, S., Acelajado, M. C., Bakris, G. L., Berlowitz, D. R., Cífková, R., Dominiczak, A. F., et al. (2018). Hypertension. Nat. Rev. Dis. Primers 4:18014. doi: 10.1038/nrdp.2018.14

PubMed Abstract | CrossRef Full Text | Google Scholar

Ottesen, J. T., Olufsen, M. S., and Larsen, J. K. (2004). Applied Mathematical Models in Human Physiology. Philadelphia, PA: SIAM.

Google Scholar

Paeme, S., Moorhead, K. T., Chase, J. G., Lambermont, B., Kolh, P., D’orio, V., et al. (2011). Mathematical multi-scale model of the cardiovascular system including mitral valve dynamics. Application to ischemic mitral insufficiency. Biomed. Eng. Online 10:86. doi: 10.1186/1475-925X-10-86

PubMed Abstract | CrossRef Full Text | Google Scholar

Pagani, E. D., Alousi, A. A., Grant, A. M., Older, T. M., Dziuban, S. W. Jr., and Allen, P. D. (1988). Changes in myofibrillar content and Mg-ATPase activity in ventricular tissues from patients with heart failure caused by coronary artery disease, cardiomyopathy, or mitral valve insufficiency. Circ. Res. 63, 380–385. doi: 10.1161/01.res.63.2.380

CrossRef Full Text | Google Scholar

Panayidou, K., Gsteiger, S., Egger, M., Kilcher, G., Carreras, M., Efthimiou, O., et al. (2016). GetReal in mathematical modelling: a review of studies predicting drug effectiveness in the real world. Res. Synth. Methods 7, 264–277. doi: 10.1002/jrsm.1202

PubMed Abstract | CrossRef Full Text | Google Scholar

Parrinello, G., Paterna, S., Torres, D., Di Pasquale, P., Mezzero, M., La Rocca, G., et al. (2009). One-year renal and cardiac effects of bisoprolol versus losartan in recently diagnosed hypertensive patients: a randomized, double-blind study. Clin. Drug. Investig. 29, 591–600. doi: 10.2165/11315250-000000000-00000

PubMed Abstract | CrossRef Full Text | Google Scholar

Pralong, W. F., Hunyady, L., Várnai, P., Wollheim, C. B., and Spät, A. (1992). Pyridine nucleotide redox state parallels production of aldosterone in potassium-stimulated adrenal glomerulosa cells. Proc. Natl. Acad. Sci. U.S.A. 89, 132–136. doi: 10.1073/pnas.89.1.132

PubMed Abstract | CrossRef Full Text | Google Scholar

Proshin, A. P., and Solodyannikov, Y. V. (2006). Mathematical modeling of blood circulation system and its practical application. Autom. Remote Control 67, 329–341. doi: 10.1134/S000511790602010X

CrossRef Full Text | Google Scholar

Reid, I. A. (1996). Angiotensin II and baroreflex control of heart rate. Physiology 11, 270–274. doi: 10.1152/physiologyonline.1996.11.6.270

CrossRef Full Text | Google Scholar

Riegger, G. A., Liebau, G., and Kochsiek, K. (1982). Antidiuretic hormone in congestive heart failure. Am. J. Med. 72, 49–52. doi: 10.1016/0002-9343(82)90576-9

CrossRef Full Text | Google Scholar

Rosalina, T. T., Bouwman, R. A., van Sambeek, M. R. H. M., van de Vosse, F. N., and Bovendeerd, P. H. M. (2019). A mathematical model to investigate the effects of intravenous fluid administration and fluid loss. J. Biomech. 88, 4–11. doi: 10.1016/j.jbiomech.2019.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Rullmann, J. A. C., Struemper, H., Defranoux, N. A., Ramanujan, S., Meeuwisse, C. M. L., and van Elsas, A. (2005). Systems biology for battling rheumatoid arthritis: application of the Entelos PhysioLab platform. Syst. Biol. (Stevenage) 152, 256–262. doi: 10.1049/ip-syb:20050053

CrossRef Full Text | Google Scholar

Runarsson, T. P., and Yao, X. (2000). Stochastic ranking for constrained evolutionary optimization. IEEE Trans. Evol. Comput. 4, 284–294. doi: 10.1109/4235.873238

CrossRef Full Text | Google Scholar

Schmitt, F., Martinez, F., Brillet, G., Nguyen-Khoa, T., Brouard, R., Sissmann, J., et al. (1998). Acute renal effects of AT1-receptor blockade after exogenous angiotensin II infusion in healthy subjects. J. Cardiovasc. Pharmacol. 31, 314–321. doi: 10.1097/00005344-199802000-00019

PubMed Abstract | CrossRef Full Text | Google Scholar

Seidlerová, J., Staessen, J. A., Maillard, M., Nawrot, T., Zhang, H., Bochud, M., et al. (2006). Association between arterial properties and renal sodium handling in a general population. Hypertension 48, 609–615. doi: 10.1161/01.HYP.0000240516.60040.ba

CrossRef Full Text | Google Scholar

Škrtić, M., Lytvyn, Y., Yang, G. K., Yip, P., Lai, V., Silverman, M., et al. (2015). Glomerular haemodynamic profile of patients with Type 1 diabetes compared with healthy control subjects. Diabet. Med. 32, 972–979. doi: 10.1111/dme.12717

PubMed Abstract | CrossRef Full Text | Google Scholar

Stanton, A., Jensen, C., Nussberger, J., and O’Brien, E. (2003). Blood pressure lowering in essential hypertension with an oral renin inhibitor, aliskiren. Hypertension 42, 1137–1143. doi: 10.1161/01.HYP.0000101688.17370.87

CrossRef Full Text | Google Scholar

Stefanadis, C. I. (2010). Complex interrelationships between heart and kidneys: establishing the role of cardiorenal syndrome. Hellenic J. Cardiol. 51, 87–88.

Google Scholar

Stéphanou, A., and Volpert, V. (2016). Hybrid modelling in biology: a classification review. Math. Model. Nat. Phenom. 11, 37–48. doi: 10.1051/mmnp/201611103

CrossRef Full Text | Google Scholar

Thomas, S. R., Abdulhay, E., Baconnier, P., Fontecave, J., Françoise, J. P., Guillaud, F., et al. (2007). SAPHIR – a multi-scale, multi-resolution modeling environment targeting blood pressure regulation and fluid homeostasis. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2007, 6649–6652. doi: 10.1109/IEMBS.2007.4353884

PubMed Abstract | CrossRef Full Text | Google Scholar

Uttamsingh, R. J., Leaning, M. S., Bushman, J. A., Carson, E. R., and Finkelstein, L. (1985). Mathematical model of the human renal system. Med. Biol. Eng. Comput. 23, 525–535. doi: 10.1007/BF02455306

PubMed Abstract | CrossRef Full Text | Google Scholar

West, G. B., Brown, J. H., and Enquist, B. J. (1997). A general model for the origin of allometric scaling laws in biology. Science 276, 122–126. doi: 10.1126/science.276.5309.122

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiśniowska, B., Tylutki, Z., and Polak, S. (2017). Humans vary, so cardiac models should account for that too! Front. Physiol. 8:700. doi: 10.3389/fphys.2017.00700

PubMed Abstract | CrossRef Full Text | Google Scholar

Wright, S. P., Moayedi, Y., Foroutan, F., Agarwal, S., Paradero, G., Alba, A. C., et al. (2017). Diastolic pressure difference to classify pulmonary hypertension in the assessment of heart transplant candidates. Circ. Heart Fail. 10:e004077. doi: 10.1161/CIRCHEARTFAILURE.117.004077

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Chen, J., Cheong, B., Pednekar, A., and Muthupillai, R. (2019). High frame rate cardiac cine MRI for the evaluation of diastolic function and its direct correlation with echocardiography. J. Magn. Reson. Imaging 50, 1571–1582. doi: 10.1002/jmri.26791

PubMed Abstract | CrossRef Full Text | Google Scholar