We utilize the comprehensive computational model of human cardiovascular and renal systems presented by Kutumova et al. (2021) and available in the BioModels database (Malik-Sheriff et al., 2020) with ID MODEL2202160001 ¹ .

A module is defined as part of a complete mathematical model with inputs and outputs to communicate with other modules:

• inputs receive variables calculated outside of the current module;

The original model (Kutumova et al., 2021) has two modular representations:

• version with two modules characterizes the interaction of the cardiovascular and renal systems;

Since the model describes the dynamics of the circulatory system, its behavior is related to the cardiac cycle. In this regard, we say that the model is in an equilibrium state if all of its variable values either do not change (e.g., systolic/diastolic blood pressure) or oscillate steadily with an amplitude equal to the length of the cardiac cycle (e.g., systemic arterial pressure).

The agent-based approach is used to model complex systems comprising a large number of subsystems, in particular, when it is difficult or impossible to formalize the system behavior at the global level. Within the framework of this approach, the system is considered as an aggregation of multiple parts—agents, each of which acts separately, using a set of rules, and can interact with other agents (Soheilypour and Mofrad, 2018; An et al., 2009; Macal and North, 2010). In our case, the agent is a mathematical model combined with an appropriate numerical solver. The simulation of agents consists in solving the underlying model until the next time point. The interaction between agents is handled by a separate scheduler and consists in the exchange of variable values between models.

In the case of the model of the cardiovascular and renal systems, the agent-based approach is utilized, since different submodels have vast differences in time scales: cardiovascular processes take fractions of a second, while the renal system contains long-term processes lasting minutes, which makes the whole model very stiff. Thus, we divide the model into two agents which are simulated with different solvers and exchange variable values during the simulation. To simulate each agent (i.e., to solve the Cauchy problem), we used a version of the CVODE solver (Hindmarsh et al., 2005) ported to Java and adapted to the BioUML application programming interface.

The Frank-Starling relationship governs normal ventricular function and ensures that the volume the heart ejects in systole (stroke volume, SV) equals the volume it receives in venous return (Costanzo, 2014). That is, within physiologic limits, the heart pumps all the blood that returns to it by the way of the veins (Hall, 2011). As venous return increases, ventricular end-diastolic volume (EDV) also increases, and due to the cardiac length-tension relationship, SV increases accordingly (Figure 2A). Note that changes in heart contractility (or inotropism) affect the Frank-Starling curves (Hillegass, 2011; Costanzo, 2014). Positive inotropic agents produce an increase in SV for a given EDV and, as a result, an increase in ejection fraction (EF), calculated as E F = S V E D V ∙ 100 %

Negative inotropic agents have the opposite effect (Costanzo, 2014). In experienced endurance athletes, the myocardial contractility is higher than in untrained healthy individuals (Schattke et al., 2014; Tanaka et al., 1986), whereas heart failure is often caused by impaired contractility (Braunwald et al., 2001). In patients with heart failure, cardiac output and external ventricular performance at rest are generally within normal limits, but an elevation of LV EDV (or preload) is associated with increases in the pulmonary capillary pressure, contributing to the dyspnea, while an elevation of right ventricular preload raises systemic venous pressure and contributes to the development of edema. LV failure becomes fatal when the SV-EDV curve is depressed to the point at which cardiac performance fails to satisfy the requirements of the body even at rest (Figure 2A).

Changes in hemodynamic parameters that occur during exercise are complex. Hyperventilation, the pumping action of the muscles, and the venoconstriction normally augment venous return. Simultaneously, the increased traffic of adrenergic nerve impulses to the myocardium, the increased concentration of circulating catecholamines, and the tachycardia augment the contractile state of the myocardium, elevating the SV without changing (or even with a reduction) of EDV (Braunwald et al., 2001). At low to moderate exercise intensity, the Frank-Starling mechanism is thought to be mainly responsible for the observed increase in SV, while submaximal and vigorous exercise involves the combination of myocardial contractility, tachycardia and the Frank-Starling mechanism (Warburton et al., 2002; Vieira et al., 2016). In heart failure, the improvement of contractility during exercise is attenuated or even prevented by norepinephrine depletion and downregulation of myocardial β-receptors, so the SV-EDV curve practically does not shift (Braunwald et al., 2001).

In our model, in accordance with the work by Proshin and Solodyannikov (2006), the maintenance of the Frank-Starling law is provided by the formula: S V = K ∙ S V max ∙ δ E D V where S V max denotes the theoretical maximum SV depending on the body weight, e.g., 200 ml for a 70 kg person, which is close to the SV values in top endurance athletes (Cumming, 1975; Saltin, 1969); δ E D V is the bell-shaped function of EDV; and K corresponds to the ventricular inotropic coefficient calculated from the ventricular inotropic state K 0 and the function σ , which characterizes sympathetic inotropic sensitivity and adaptive capacity of the myocardium depending on the neurohumoral factor H: K = K 0 + 0.25 ∙ σ H , K 0 .

K 0 is a constant parameter that takes values from 0.2 to 0.8. The exact form of the functions δ and σ is determined in the supplementary file of our base article (Kutumova et al., 2021), so for simplicity we don’t detail them here. Figure 2B shows an example of the relationship curves between LV SV and EDV simulated by the model. In this example, we considered the equilibrium model parametrization reproducing the physiological characteristics of a normal adult with K L 0 = 0.55 . To bring the model out of equilibrium and obtain Frank-Starling curves, we changed the R O 2 parameter, which is responsible for the body’s oxygen demand, from the normal value of 4.2 ml/s (≈250 ml/min) (Nathan and Singer, 1999) to 16.7 ml/s (≈1,000 ml/min), i.e., the level of moderate exercise (Thomson, 1971; Julius et al., 1967). Thus, we obtained the SV-EDV curve before the model enters a new equilibrium state at K L 0 = 0.55 . As this parameter was increased, the curve shifted upward and to the left, while decrease in K L 0 resulted in a shift downward and to the right. Critically low contractility values lead in the model to a persistent increase in pulmonary venous pressure, which after some model time results in a sharp drop in systolic pressure and SV (Figure 2C). We interpret this state as cardiogenic shock (Hollenberg et al., 1999; Handler, 1985), since hemodynamic criteria of this clinical condition are sustained hypotension (systolic blood pressure < 90 mmHg for at least 30 min) and a reduced cardiac index (<2.2 L/min per m²) in the presence of elevated pulmonary capillary occlusion pressure (>15 mmHg) (Hollenberg et al., 1999), which serves as an estimate of pulmonary venous pressure (Grignola, 2011; Chaliki et al., 2002). An interactive implementation of this example in BioUML is available online (see the Availability section).

We performed drug activity modeling as follows. First, we defined target variables in the model for all antihypertensive agents, and then we multiplied each target variable by the specific function F that provides the pharmacodynamic effect of therapy: F = 1 ± ∑ D E D ∙ d F (1)

In this equation, E D denotes the magnitude of target stimulation (the plus sign) or target inhibition (the minus sign) as a function of drug dose D. The summation is made for all drugs belonging to the same antihypertensive class and for all simulated dosages. The drug indicator d F is equal to 1 or 0 depending on whether a patient is being treated with this drug or not. All pharmacodynamic effects are initialized in a special module “Pharmacodynamics”. If only one indicator is equal to 1, then the patient receives monotherapy, if not, then combination therapy. The targets for the different antihypertensive agents incorporated in the model are shown in Figure 1 and listed in Table 1. The values of all sums ∑ D E D ∙ d F , denoted in the pharmacodynamics module as A C E i (angiotensin-converting enzyme inhibition), A R B (angiotensin receptor blocking effect), D R I (direct renin inhibition), C C B e a , C C B a a (calcium channel blocking effects), etc. (Figure 1), are passed from this module to the top level of the model and further to the modules where they are required. Below we provide a description of the pharmacodynamic functions for all targets.

Aliskiren is a low molecular weight renin inhibitor of non-peptide structure (Gradman et al., 2005). The antihypertensive effect, 1 – D R I , is implemented in the module “Angiotensin” by reducing plasma renin activity. The target inhibition is determined using the E _max -model (Derendorf and Meibohm, 1999): E D = E max ∙ D D + E D 50 where E max = 0.99 is the maximum drug effect, and E D 50 = 20 mg is the concentration of the medicine that gives half of the maximal effect (Hallow et al., 2014). In the model, we consider two dosages of aliskiren, 150 and 300 mg. If we denote the corresponding indicators of the drug as A l i s k i r e n 150 and A l i s k i r e n 300 , then the value of D R I is calculated by the formula: D R I = E 150 ∙ A l i s k i r e n 150 + E 300 ∙ A l i s k i r e n 300

Amlodipine causes dilation (and, as a result, a decrease in resistance) of afferent and efferent arterioles (Hayashi et al., 2007), as well as interlobar, arcuate, and interlobular arteries (Hallow et al., 2014). The corresponding influences, 1 – C C B a a , 1 – C C B e a , and 1 – C C B p r e g l o m , appear in the module “Glomerular filtration”. In addition, amlodipine acts on systemic vascular resistance (Hallow et al., 2014), which primarily depends on the conductivity of systemic microvessels (arterioles, capillaries and venules). This effect, 1 – C C B s y s , is used in the module “Systemic arteries'', where it is applied to the blood flow from the arterial to the venous bed of the systemic circulation.

Bisoprolol is a β1-adrenoceptor antagonist with no partial agonist (intrinsic sympathomimetic) activity or membrane stabilizing (local anaesthetic) activity (Lancaster and Sorkin, 1988). Since β-blockers suppress renin secretion in the kidneys and reduce its level and activity (Laragh and Sealey, 2011), we added the corresponding effect, 1 – B b l o c k e r _ r s , into the module “Renin''. These drugs also block the access of catecholamines (norepinephrine, epinephrine) to their receptors so that the heart rate and blood pressure are reduced (Frishman, 2003). Therefore, we introduced the influence of 1 – B b l o c k e r on the activity of stress receptors in the module “Neurohumoral control”.

Enalapril is an angiotensin-converting enzyme inhibitor (Ferguson et al., 1982) which acts in the module “Angiotensin''. The action of the drug is modeled by multiplying the rate of conversion of angiotensin I to angiotensin II by the function 1 – A C E i (Hallow et al., 2014).

Hydrochlorothiazide (HCTZ) belongs to the benzothiadiazine class, referred to simply as thiazide diuretics (Carter et al., 2004). Thiazide-induced reduction of arterial pressure includes differentiation into acute and chronic phases (Rapoport and Soleimani, 2019). The acute reduction correlates with diuresis and a decrease in plasma volume associated with inhibition of sodium reabsorption in the distal tubules (Duarte and Cooper-DeHoff, 2010). This phenomenon is described by the drug effect of 1 – D i u r e t i c I n h i b i t i o n on the normal value of fractional distal tubule sodium reabsorption (Hallow et al., 2014) in the module “Sodium”.

A possible mechanism of the chronic reduction includes vascular dilation (Rapoport and Soleimani, 2019; Duarte and Cooper-DeHoff, 2010). Thus, by analogy with the model by Hallow et al. (2014), we introduced into the module “Glomerular filtration'' the effects of 1 – D i u r e t i c a a , 1 – D i u r e t i c e a , and 1 – D i u r e t i c p r e g l o m on the resistances of afferent arterioles, efferent arterioles, and interlobar/arcuate/interlobular arteries, respectively. Considering the similar effect on systemic microvessels, we took into account the following facts:

• extracellular fluid (as well as plasma volume) is initially reduced by thiazides and then almost fully recovers within 4–6 weeks of continuous treatment (Duarte and Cooper-DeHoff, 2010; Rapoport and Soleimani, 2019; Leth, 1970);

We found that the recovery of these values can be reproduced in the module “Systemic arteries” with a growing influence of 1 / 1 – D i u r e t i c s y s on the conductivity of systemic microvessels (the inverse of their resistance): D i u r e t i c s y s t i m e = min E max ∙ t i m e d u r a t i o n , E max where E max denotes the maximum influence and duration is the time to reach E max .

As other targets of HCTZ in the model, we considered constant levels of urea and potassium in the blood. Since HCTZ increases blood urea and decreases blood potassium (Scaglione et al., 1992; Scaglione et al., 1995; Devineni et al., 2014), we multiplied these parameters by the factors 1 + D i u r e t i c u r e a and 1 – D i u r e t i c p o t a s s i u m (the modules “Hormonal system” and “Aldosterone”, respectively).

In addition, a clinical trial by Villamil et al. (2007) shows that HCTZ monotherapy increases plasma renin activity. Therefore, following the model by Hallow et al. (2014), we took into account the direct effect of 1 + D i u r e t i c S t i m u l a t i o n on renin secretion (the module “Renin”).

Another target of the drug is based on the fact that heart rate does not change significantly with long-term treatment with HCTZ (van Brummelen et al., 1980; Shah et al., 1978; Scaglione et al., 1992; Scaglione et al., 1995). Different studies show that norepinephrine pressor response decreases or remains unchanged with thiazide diuretics (Rapoport and Soleimani, 2019). Thus, we assumed the effect of 1 – D i u r e t i c s t r e s s on the activity of stress receptors in the module “Neurohumoral control”.

Losartan is an angiotensin-receptor antagonist without agonist properties (Sica et al., 2005). Its action is modeled in the module “Angiotensin” by the function of the form 1 – A R B multiplied by the rate of angiotensin II binding to AT1 receptors (Hallow et al., 2014).

Similar to the study by Cheng et al. (2017), we define a virtual patient as a single equilibrium parameterization of the model. This parameterization corresponds to a specific state of the patient, for example, “sick” or “healthy”. Response of the patient to external stimuli (e.g., taking antihypertensive drugs) is carried out by introducing a perturbation into the model values and subsequently solving the modified system to search for a new equilibrium. We divide the model parameters (constant values of the model) into two groups:

• the personal parameters vary within normative or pathological ranges or take a fixed value in accordance with clinical measurements in a real patient;

Among the model variables, we single out observable (i.e., clinically measurable) characteristics. Such variables can be used to verify how well a virtual patient matches a real person.

The use of optimization methods in the current work was necessary in two problems:

• to calibrate the kinetic parameters of drugs based on experimental studies;

In both cases, we fitted the model parameters by minimizing the distance function, defined as the normalized sum of squared differences (Hoops et al., 2006) between several observed ( X i exp ) and simulated equilibrium ( X i ) targets (blood pressure, heart rate, etc.): f d i s t X 1 , … , X n = ∑ i = 1 n ω min ω i X i − X i exp 2 ,   ω i = X i exp ,   ω min = min ω i ,   i = 1 , … , n (2)

Additionally, to take into account the known physiological constraints Y j min ≤ Y j t ≤ Y j max , j = 1 , … , m , we considered the penalty function (Runarsson and Yao, 2000): f p e n a l t y Y 1 , … , Y m = ∑ t ( ∑ j = 1 m max ⁡ 0 , Y j min − Y j t 2 + ∑ j = 1 m max ⁡ 0 , Y j t − Y j max 2 ) (3)

To solve such optimization problems, we used a stochastic ranking evolution strategy (Runarsson and Yao, 2000) suitable for constrained global optimization.

The task of generating a virtual patient is to find an equilibrium state of the model within the given physiological ranges. Suppose we have a medical history that includes the results of laboratory diagnostics for some person (blood pressure monitoring, electrocardiography, echocardiography, blood tests, etc.). We can directly substitute those of the measured values which are constant in the model (these are body weight and blood counts). The measured values calculated in the model (such as blood pressure, heart rate, glomerular filtration rate, etc.) can be achieved by minimizing the objective function (2). If the patient’s history confirms the presence (or absence) of any cardiovascular or renal diseases, we can also consider the constraints imposed on a number of model parameters and variables that were not explicitly measured. For example, the normal range of systemic vascular resistance is 700–1,600 dyn×s/cm⁵ (Klingensmith et al., 2016), while in patients with congestive heart failure, the values are generally higher: 944–2,209 dyn×s/cm⁵ (Laskey et al., 1990). In the case of the model parameters, such constraints determine the boundaries of the search space (the lower and upper limits of the fitting parameters). In the case of variables, they define the feasible region used to calculate the penalty function (3). A complete list of the ranges of parameters and variables that we consider in the current work is provided in Supplementary Tables S1, S2 and justified in our basic study (Kutumova et al., 2021).

Note that the agent-based model converges to the equilibrium state much more slowly than each of the agents independently. Therefore, the optimization of the entire model takes more time. However, the mean arterial pressure, cardiac output and hematocrit values sent by the cardiovascular submodel to the renal agent (Figure 3) are clinically measurable and can be found in the patient’s history (if not, we can randomly select them from physiological ranges). By fixing these values, we can estimate the parameters of only the renal system, pass the necessary equilibrium values to the cardiovascular agent, and then optimize its parameters by minimizing the distances between the simulated and fixed equilibriums of mean arterial pressure and cardiac output at a constant value of hematocrit. Figure 3 demonstrates the described algorithm.

In the current work, we generate virtual patients with uncomplicated hypertension. Analyzing the model, we found that some equilibrium parameterizations under perturbations of physiological parameters result in persistent increase in pulmonary venous pressure, which occurs due to the accumulation of blood in the pulmonary veins. This condition leads to pulmonary edema (Burkhoff and Tyberg, 1993; Alwi, 2010), which may be the result of a hypertensive crisis (Varounis et al., 2017). Because of the complexity of the model, the regions in the parameter space that induce such a state are not obvious and represent a separate issue for study, which is beyond the scope of this work. Here, we restrict ourselves to excluding invalid patients from the analysis. An increase in pulmonary venous blood volume can be detected with an increase in total blood volume that follows an increase in sodium intake (denoted as Φ s o d i n in the model). Therefore, we checked each generated virtual patient by simulating the following sodium loading test. We examined the instantaneous change in Φ s o d i n from a patient value in the range of 0.0280–0.2088 mEq/min to the elevated value of 0.243 mEq/min (He et al., 2001), and then observed the subsequent behavior of the model variables. Acceptable dynamics was characterized by convergence to a new equilibrium state with an increased systolic blood pressure (SBP) and a slightly increased pulmonary venous pressure. We considered 25 mmHg to be the maximum possible elevation of SBP based on the experimental data (He et al., 2001) reporting a mean change in this parameter of approximately 17 mmHg between low and high sodium intake in hypertensive patients.

A virtual population is a collection of unique virtual patients. To generate a population with varying values of systolic blood pressure, diastolic blood pressure (DBP), heart rate (HR), body weight, and body mass index (BMI), we used a function producing the random numbers with the required mean and standard deviation (SD).

We analyzed the baseline characteristics of hypertensive patients enrolled in different clinical trials of aliskiren, amlodipine, bisoprolol, enalapril, hydrochlorothiazide, and losartan (Abate et al., 1998; Agabiti Rosei et al., 2005; Brown et al., 2011; Corea et al., 1996; Cushman et al., 2012; Derosa et al., 2014; Fagher et al., 1990; Fermé et al., 1990; Gradman et al., 2005; Koh et al., 2004; Kwakernaak et al., 2017; Lee et al., 2012; Lithell et al., 1987; Mallion, 2007; Martina et al., 1999; Narkiewicz, 2007; Nedogoda et al., 2013; Oparil et al., 1996; Oparil et al., 2008; Pareek et al., 2016; Paterna et al., 2007; Pool et al., 2007; Porthan et al., 2009; Puig et al., 2007; Sun et al., 2016; Tham et al., 1993; Weber et al., 1995; Wing et al., 2003). Based on these characteristics, we selected suitable mean and SD values for the necessary quantities. A summary table describing the clinical trial populations is given in the Supplementary Material (Supplementary Table S3). Table 2 shows the settings of SBP, DBP, HR, BMI, weight, height, and gender that we used in our work. Note that for other parameters and variables of the model, we also took into account physiological ranges typical for a patient with uncomplicated arterial hypertension (Supplementary Tables S1, S2). Thus, as the inclusion criteria, we considered essential hypertension (SBP ≥ 140 mmHg and/or DBP ≥ 90 mmHg; SBP > 130 mmHg; DBP > 80 mmHg), BMI > 22 kg/m², and average height (160–180 cm). The exclusion criteria were severe hypertension (SBP > 179.5 mmHg or DBP > 109.5 mmHg), cardiovascular or renal complications, abnormal heart rate (HR < 60 or HR > 90 beats/min), and severe obesity (BMI > 36 kg/m²).

The algorithm for obtaining virtual patients was as follows:

1. Randomly generate SBP, DBP, HR, BMI, and weight values with the basic settings from Table 2 so that these values do not fall into the exclusion criteria.

We performed the model calibration in the following way. For the pharmacodynamics of aliskiren, we used the kinetic parameters from the study by Hallow et al. (2014). Since aliskiren does not affect HR (Stanton et al. 2003) and glomerular filtration rate (GFR) (Kwakernaak et al., 2017; Bokuda et al., 2018; Siddiqi et al., 2011), we reproduced the similar behavior of these variables in the model. For this purpose, we generated a test population of 100 virtual hypertensive patients with varying values of SBP, DBP, HR, body weight, and BMI (as described in the current section above). Then, for each patient in the population, we estimated the dimensionless parameters responsible for the effect of angiotensin II on the cardiac muscle through baroreceptors ( s l b a r o ) and stress receptors ( s l s t r e s s ), as well as the effect of angiotensin II on the glomerular filtration coefficient ( s l K F G ) (Kutumova et al., 2021). We performed all estimations so that the equilibrium values of HR and GFR before treatment with aliskiren coincided with the simulated equilibrium values after treatment. Next, we found the mean values of the estimated parameters s l b a r o , s l s t r e s s and s l K F G in the population: s l b a r o m e a n = 0.0499 ,   s l s t r e s s m e a n = 0.2133 ,   s l K F G m e a n = 0.2015 , and fixed them in the model. This change led to a violation of the equilibrium states of the virtual patients. Therefore, at the next step, we found new equilibrium states. After that, we fitted the kinetic parameters of all other drugs for each patient in the adjusted population so that the simulated changes in the target physiological variables (given in the last column of Table 1) were equal to the mean changes in the corresponding experimental characteristics. As the final values of the therapy parameters [column “E(D)-value” in Table 1], we took the mean values for the population.

BioUML (https://ict.biouml.org/; https://sirius-web.org/bioumlweb/) is a Java-based integrated environment for systems biology (Kolpakov et al., 2022). It supports a wide range of biological formats and such tools as visual modeling, simulation, parameter estimation and a number of other numerical methods. The software features used in this work are the following:

• a web-version (for collaboration and public presentation of data) and a standalone version (for independent work) of the program;

One of the techniques for generating virtual patients is described by Allen et al. (2016). The authors propose to create a large number (hundreds of thousands) of “plausible patients”, defined as model parameterizations within biologically plausible ranges, and then select from them a virtual population corresponding to the desired empirical distribution. To find a plausible patient, they take model inputs within the predefined plausible bounds and optimize this choice until the required outputs also fall within plausible ranges. However, the authors note that their approach is not applicable to models that are slowly simulated (for example, with dynamics across multiple time-scales, as is the case of our model) due to high computational cost to get a large plausible population. An easier way to create a virtual population is discussed in works by Hallow et al. (2014), Hallow et al. (2020), where inputs are sampled from predefined ranges and considered as valid virtual patients only if they produce outputs with physiologically reasonable values. This approach may be relevant with a relatively small number of inputs and outputs [the authors consider 11 inputs/13 outputs in one study (Hallow et al., 2014) and 14 inputs/3 outputs in another study (Hallow et al., 2020)]. However, with a large number of inputs/outputs (in our case, 57 inputs given in Supplementary Table S1 and 54 outputs given in Supplementary Table S2), it leads to enumeration of a huge number of samples with a high rejection rate. Therefore, we apply our own method for generating a virtual population. Like Allen et al., we use optimization methods with a penalty function to create virtual patients, which ensures that outputs are within physiologically acceptable ranges. In addition, for each patient, we minimize the objective function of the distances between the simulated and normally distributed random values of SBP, DBP, and HR, which directly gives us a population with predefined empirical distributions of these quantities (for more details see the Materials and methods section). To test the model, we generated a population of 250 virtual hypertensive patients. The baseline characteristics of this population are shown in Table 3.

Using the pharmacodynamic models determined above, we simulated 4-week antihypertensive treatment of the generated virtual patients. To check the resulting changes in the model variables, we compiled a table describing the physiological response to therapy with aliskiren, amlodipine, bisoprolol, enalapril, HCTZ, and losartan (Supplementary Table S4). This table combines information from previously published clinical studies and includes post-treatment changes (increase, decrease, or no change) in plasma renin activity, plasma concentrations of renin, angiotensin I, angiotensin II, aldosterone, serum sodium concentration, heart rate, cardiac output, glomerular filtration rate, filtration fraction, renal blood flow, renal vascular resistance, extracellular fluid volume, and afferent arteriolar resistance. An increase (or decrease) in the long-term post-treatment dynamics of a variable was considered achieved if the corresponding change in value exceeded 5%. At the same time, if the change was less than 10%, we regarded it as insignificant. Thus, the range of changes in dynamics from 5 to 10% was taken as a transitional range due to the variability of possible changes in population variables. Virtual patients with the appropriate response (186 subjects out of 250, or 74%, “Patients included in analysis” in Table 3) were selected for further analysis. Figure 4 shows the distribution of the main parameters in this population. Distribution plots for more parameters are available in the Jupyter document on the web-version of BioUML (see the Availability section).

The results of the comparison of the simulated changes in systolic and diastolic blood pressure with experimental measurements are shown in Figure 5. As can be seen from these data, the model accurately reproduces the reduction in blood pressure for all drugs. Numerical data related to Figure 5 and the Jupyter histogram visualization can be found in the web-version of BioUML (see the Availability section).

In our basic study (Kutumova et al., 2021), we showed the model’s ability to simulate both healthy subjects and patients with cardiovascular pathologies, including systemic arterial hypertension, heart failure, pulmonary hypertension, etc. In the current work, we introduce the algorithm for generating virtual populations. The integration of these capabilities allows the creation of virtual populations with various combinations of cardiovascular diseases and the design of in silico experiments depending on the purpose of the investigation. A common problem in drug research is comparing the effects of different medicines and doses in different categories of patients. Below we provide two examples of such a comparison for mono- and combination therapies performed in the virtual population described in the previous section. All graphs given in the examples were implemented as Jupyter documents in BioUML (see the Availability section for details).

The main limitation of the work is related to the generation of a virtual population. We consider a large number of constraints imposed on the parameters and variables of the model (Supplementary Tables S1, S2). But this is not enough to cut off all virtual patients with unrealistic dynamics. As a solution to the constraint optimization problem, we can obtain a parameterization of the model that demonstrates an excessive increase in systolic blood pressure (more than 25 mmHg) with an increase in sodium intake to a value of 0.243 mEq/min, or parameterization that leads to incorrect model dynamics when simulating antihypertensive therapy (for example, an increase in heart rate during treatment with bisoprolol). In the current work, we removed such parameterizations with additional checks. However, their isolation at the optimization stage could significantly speed up the process of generating virtual patients. One of the ways to solve this problem is to look for correlations between the domain (and/or range) of the model and the resulting dynamics of control variables in response to external stimuli (a sharp increase in sodium intake, activation of one of the antihypertensive treatment regimens, etc.). Adding appropriate constraints (determined from such correlations) to the optimization problem will result in valid virtual patients. However, the values of the correlation parameters must be close to unity, otherwise we risk losing some of the valid solutions. The second possible way to solve this problem is to introduce simulation experiments (test cases) into the optimization problem. In this case, each potential solution at all iterations of the optimization algorithm will be checked not only for the fulfillment of the established parametric constraints, but also for completion of test cases. This approach avoids the loss of solutions, as in the case of insufficiently high values of correlation parameters, but requires a change in the system implementation of optimization algorithms in BioUML.

Another limitation of the study is the complexity of the considered model, which is directly related to the cost of resources and time and is reflected in the speed of generation of virtual patients. However, the constrained optimization algorithm used to obtain a single patient is easily parallelized, as is the process of creating a large population. Therefore, with the development of computing technologies and the emergence of more powerful servers with multi-core processors, this problem will become less acute. Low-dimensional models certainly have a number of benefits besides saving time and enhancing productivity. In particular, they are easier to understand and interpret. And in some cases, they reproduce experimental data and make reasonable predictions just as well as more complex models (Kutumova et al., 2013). Therefore, simple mechanistic models are still useful and relevant for solving applied problems related to the study of pathological conditions of human physiology. However, the main disadvantage of simple models is the loss of biological information. Thence, we believe that over time, physiological models will continue to evolve towards greater complexity and closer to reality.

We consider each parameterization of the model in an equilibrium state and within physiological ranges as a virtual patient. To relate this concept to a real patient, we can adjust the model parameters so that the equilibrium dynamics reproduces laboratory measurements of desired physiological characteristics (obtained from blood pressure monitoring, electrocardiography, echocardiography, blood tests, etc.). However, these data give only a small part of the model quantities. To control the unknown parameters and variables, we can move from a single virtual patient to a virtual population and take into account significant variation in model values. Treatment simulation of such a population allows it to be divided into groups with a similar response to drugs. Thus, we can find out which patient characteristics contribute to the effectiveness (or failure) of antihypertensive therapy. Linking the model to real patients is the main direction of our future work.