The original model (Karatza and Karalis, 2020) was reproduced in the BioUML ⁴ software (Kolpakov et al., 2022; Kolpakov et al., 2019) using SBML format for technical representation (Keating et al., 2020) and SBGN format for visualization (Gambardella et al., 2009) (Figure 2).

The model includes four compartments: stomach, small intestine, central, and peripheral compartments. These compartments contain two substances, losartan and E-3174, which are interconnected through first-order reaction equations (numbers 1 through 7 in Figure 2). A detailed description of all equations and variables is given in Supplementary Tables S1, S2, respectively.

The model also includes discrete events, initial assignments, and algebraic and differential equations used to solve the model. The list of mathematical tools and variables can be found in Supplementary Tables S3, S4, respectively.

The pharmacokinetics of losartan were assessed by calculating the following parameters: area under the concentration-time curve (AUC) for losartan (AUC_losartan) and its metabolite (AUC_E-3174), AUC_losartan to AUC_E-3174 ratio (AUC_ratio), maximum plasma concentration (C_max), time to reach C_max (t_max), and terminal elimination half-life (t_1/2) (Table 3).

C_max and t_max were derived from the model simulation results. AUC_losartan and AUC_E-3174 were calculated using the differential equations: d AUC losartan d t i m e = C p , d AUC E − 3174 d t i m e = C m , where C_p and C_m denote concentrations of losartan and E-3174, respectively.

Since the model was implemented using delay differential equations in the conversion of losartan to E-3174, the concentration of E-3174 in the first numerical steps of the simulation is zero. Consequently, the AUC_E-3174 value is also zero. Therefore, a piecewise function was used to calculate AUC_ratio: AUC ratio = 0 , if  C m = 0 , AUC losartan AUC E − 3174 , otherwise

The t_1/2 value was calculated by linear regression analysis from the terminal linear part of the plasma concentration versus time semilogarithmic plot according to the following equation: t 1 / 2 = t 2 − t 1 × ln 2 ln C 1 C 2 where C₁ and C₂ are the concentrations of the corresponding substance at times t₁ and t₂, respectively, according to the linear regression line (Toutain and Bousquet-Mélou, 2004).

To test the losartan metabolism model in a population of virtual patients with different CYP2C9 genotypes, we used the computational model of the human cardiovascular and renal systems that, in particular, simulates the pharmacological effects of losartan on cardiovascular and renal parameters (Kutumova et al., 2022; Kutumova et al., 2021). The model is discrete-continuous and consists of a system of ordinary differential equations with several discrete events corresponding to instantaneous changes in the modeled dynamics (e.g., the transition from systole to diastole). This model is available in the BioModels database (Malik-Sheriff et al., 2020) with ID MODEL2202160001 ⁵ .

A virtual patient is a single equilibrium parameterization of the model within physiological limits. In the current work, we used 100 virtual hypertensive patients from our previous study on antihypertensive therapy modeling (Kutumova et al., 2022).

To simulate treatment of virtual patients with losartan, we used the following equation from the cardiorenal model (Kutumova et al., 2022), which reduces the rate of angiotensin II binding to AT1-receptors: A R B = k b l o c k × L o s a r t a n t r e a t m e n t , where ARB (angiotensin receptor blocking effect) is the AT1-receptor blocking activity of the drug, Losartan_treatment is a discrete parameter that can take values of 0 (no treatment) or 1 (treatment course), and k _block is the total AT1-receptor blocking activity of losartan across CYP2C9 genotypes.

The value of the parameter k _block is fed into the cardiorenal model from the losartan metabolism model and depends on the CYP2C9 allelic variant. Parameter values for the different CYP2C9 genotypes were estimated using the following considerations. Losartan 25 mg orally once daily has not been shown to produce clinically significant reductions in blood pressure compared with placebo (Gradman et al., 1995). A similar effect in the cardiorenal model is given by k _block = 0.1. In addition, for daily doses of 50 and 100 mg losartan, k _block values were estimated to be 0.886 and 0.954 (Kutumova et al., 2022). Using these doses in the pharmacokinetic model with parameters fitted to the CYP2C9∗1/CYP2C9∗1 genotype (wild-type), we calculated the corresponding AUC_E-3174 values. Thus, we received three relationship points between k _block values and corresponding AUC_E-3174 values (Supplementary Table S5).

To describe this dependence, we used the following E-max model (Kirby et al., 2011): k b l o c k = E max × AU C E − 3174 α E D 50 α + AU C E − 3174 α

The coefficients of the fitted E-max model are presented in Supplementary Table S6. Finally, using this equation, we extended the relationship between k _block and AUC_E-3174 values to any values of AUC_E-3174, i.e., for any CYP2C9 allelic variant (Table 4) (Supplementary Figure S1).

To simulate the models, we used a version of the CVODE solver (Hindmarsh et al., 2004) ported to Java and adapted to the BioUML software interface.

CVODE solves initial-value problems for systems of ordinary differential equations in real N-space. It is used to solve both stiff and non-stiff systems by variable-order and variable-step multistep methods. CVODE contains two groups of multistep formulas that are suitable for different systems. For non-stiff problems, CVODE includes the Adams-Moulton formulas, and for stiff problems, the backward differentiation formulas. The coefficients for these methods are determined based on the method type, its order, the recent history of step sizes, and the normalization α_n,0 = −1.

To solve the inverse problem of identifying model parameters based on experimentally measured variables, nonlinear optimization methods were used.

This problem was solved by finding the minimum of the objective function, which was calculated by summing the squares of the distances from the experimental points to the points predicted by the model. If the problem required additional constraints imposed on some variables of the model, a penalty function was used to assess their feasibility.

BioUML contains the following optimization methods:

1. Stochastic Ranking Evolutionary Strategy (SRES) (Runarsson et al., 2000),

Three of these methods (SRES, MOPSO and MOCell) were most suitable for our task. We compared the optimization results of these methods for our problem and found that they differed by no more than 5%. Therefore, we arbitrarily selected one of them (SRES) for further calculations.

After estimating the model parameters based on experimental data, it is important to understand how accurately these parameters have been estimated in terms of the quantity and quality of the data. This understanding is necessary for further investigation of model predictions and can be provided by analyzing the parameters for identifiability (Raue et al., 2010; Raue et al., 2009). To study the sensitivity of the objective function to changes in a fitting parameter, we exclude it from the optimization process with a fixed value that gradually increases and then decreases compared to the optimal solution. In this way, we determine the influence of this parameter on the value of the objective function (i.e., the quality of the experimental data approximation). If the shift of the parameter in any direction along the numerical axis leads to a significant increase in the objective function, then it is identifiable. If a significant increase in the objective function occurs when moving in only one direction, then the parameter is partially identifiable. Otherwise, it is impossible to determine the parameter based on the available experimental data, that is, it is unidentifiable.

To train the losartan metabolism model, we needed experimental plasma concentrations of the drug and its active metabolite in individuals with homozygous CYP2C9 genotypes (CYP2C9*1/CYP2C9*1, CYP2C9*2/CYP2C9*2, and CYP2C9*3/CYP2C9*3). We were able to find only one article with the time courses of these concentrations after a single oral dose of 50 mg losartan potassium for all the above genotypes (Yasar et al., 2002).

To digitize this data from the figures in the original study, we used Plotdigitizer ⁶ software. In analyzing the data, we decided not to use the first losartan point of the concentration-time curve from the CYP2C9*2/CYP2C9*2 plot because immediately after drug administration (time = 0), the concentration of losartan was 860 nM, which means that the absorption of the drug occurred instantly and is unlikely since it takes some time to get from the stomach to the small intestine and then into the blood. In addition, the maximum concentration values of losartan and E-3174 in the figures were replaced with the tabulated values of C_max for the corresponding genotypes from the same study to obtain more plausible results. The results of data digitizing are presented in Supplementary Table S7.

To develop and analyze the model, we used the BioUML software (version 2023.3), a Java-based integrated environment for the modeling of different biological systems (Kolpakov et al., 2022; Kolpakov et al., 2019).

The original model of losartan pharmacokinetics was validated using experimental data on the time-dependent plasma concentrations of losartan and its metabolite after a single oral dose of 100 mg losartan potassium (Figure 2) (Karatza and Karalis, 2020). At the same time, in the clinical study of individuals with different genotypes of CYP2C9, which plays a major role in the metabolism of the drug, a dosage of 50 mg was used (Yasar et al., 2002). Directly changing the drug dosage from 100 to 50 mg in the model showed a significant discrepancy between the experimental data and the model predictions (Supplementary Figure S2). Furthermore, the original model did not take into account the different CYP2C9 activities. Thus, to refine the model to incorporate CYP2C9 genetic factors, it was necessary to recalibrate its parameters.

By analyzing the list of these parameters (Supplementary Table S2), we came to the conclusion that the value of the parameter b in the sinusoidal equation (number 1 in Figure 2), which describes the periodic opening and closing of the pyloric valve of the stomach (Supplementary Table S1), has a physiological justification. It consists in the fact that the cycle of the migrating motor complex (a regular pattern of gastric motility during fasting) (Deloose et al., 2012) varies on average from 85 to 115 min (Feher, 2017), and at b = 3.95 (the value in the original model), the resulting period falls within this time interval.

In addition, the parameter k _m from reaction 3 in Figure 2 is a first-order kinetic constant characterizing the rate of conversion of losartan to E-3174, i.e., the activity of CYP2C9. Therefore, we estimated the unique value of this parameter for each of the three homozygous genotypes.

The final k _m values for each genotype are presented below: k m C Y P 2 C 9 * 1 / C Y P 2 C 9 * 1 = 1.056   h − 1 , k m C Y P 2 C 9 * 2 / C Y P 2 C 9 * 2 = 0.823   h − 1 , k m C Y P 2 C 9 * 3 / C Y P 2 C 9 * 3 = 0.025   h − 1 .

The optimized values of all other parameters are listed in Supplementary Table S8.

The time-dependent concentrations of losartan and E-3174 after a single oral dose of 50 mg losartan potassium for three CYP2C9 homozygotes predicted by the model after parameter redefinition, as well as experimentally obtained data for individuals with the same genotypes (Yasar et al., 2002), are shown in Figure 3. Because the concentration of losartan was significantly higher than the concentration of its metabolite for CYP2C9*3/CYP2C9*3 genotype, we plotted profiles of both substances separately (Figure 4). The predicted curves accurately describe the clinical data.

After validation of the model, we calculated the values of the pharmacokinetic parameters (AUC, AUC_ratio, C_max, t_1/2, and t_max) (Table 3) for the simulated profiles of losartan and E-3174 in homozygous genotypes and compared them with the corresponding values calculated from the experimental data (Table 5).

The values of all simulated pharmacokinetic parameters fell within the corresponding experimental range (Table 5). However, the t_{1/2, E-3174} value for the CYP2C9*1/CYP2C9*1 case was at the lower bound of the range. t_max values for losartan and E-3174 for different genotypes were not provided in the clinical study (Yasar et al., 2002). Therefore, we derived these values directly from the losartan and E-3174 concentration-time profiles reported in the article for homozygous CYP2C9 genotypes. Note that we were only able to derive mean values for t_max from these data, while the variation among patients was unknown because blood samples for analysis were taken at fixed times (0, 0.5, 1, 1.5, 2, 4, 6, 8, 10, 12, 24 h). In addition, the group of patients with the CYP2C9*3/CYP2C9*3 genotype included only one individual (Yasar et al., 2002) and direct calculation of SD for all pharmacokinetic parameters in this case was also not possible. For the comparative analysis of experimental and simulated values, we assumed that the unknown SD values mentioned above could be the same as for other parameters in the CYP2C9*1/CYP2C9*1 group (SD varies from 28% to 73% for different parameters, i.e., is 50% on average). Considering SD = 50% for all parameters in CYP2C9*3/CYP2C9*3 and for t_max parameters in all groups, we found that all model-predicted values, except AUC_ratio for the CYP2C9*3/CYP2C9*3 case, fell within the mean ± SD range, while the simulated value AUC_ratio = 13.5 deviated from the experimental value of 8.9 by 51%.

After optimizing the model parameters, we checked them for identifiability to ensure that the resulting solution is unique. As can be seen from Figure 5, all model parameters are identifiable.

The model was tested using the data from the same experimental study (Yasar et al., 2002), but obtained for patients with heterozygous CYP2C9 genotypes (CYP2C9*1/CYP2C9*2, CYP2C9*1/CYP2C9*3, and CYP2C9*2/CYP2C9*3).

Since the rate of conversion of losartan to E-3174 in the model was described using a first-order reaction equation, the arithmetic mean of the k _m values for the two homozygotes whose alleles were included in the heterozygote was used as k _m for this heterozygote: k m C Y P 2 C 9 * 1 / C Y P 2 C 9 * 2 = 1.056 + 0.823 2 = 0.940   h − 1 , k m C Y P 2 C 9 * 1 / C Y P 2 C 9 * 3 = 1.056 + 0.025 2 = 0.541   h − 1 , k m C Y P 2 C 9 * 2 / C Y P 2 C 9 * 3 = 0.823 + 0.025 2 = 0.424   h − 1 .

In the study by Yasar et al. (2002), there were no concentration-time curves for patients with heterozygous CYP2C9 genotypes, so it was impossible to check how well the concentration-time curves of losartan and E-3174 simulated by the model matched the experimental ones (Supplementary Figure S3). Therefore, to verify the model, we compared pharmacokinetic parameter values from the clinical study Yasar et al., (2002) and those predicted by the model (Table 6).

Because concentration-time profiles of losartan and E-3174 were not available for heterozygous CYP2C9 genotypes, we were unable to estimate t_{max, losartan} and t_{max, E-3174} values and instead used the corresponding values obtained for the CYP2C9*1/CYP2C9*1 case.

Table 6 demonstrates that the model predictions for heterozygous CYP2C9 genotypes fit the experimental data worse than for homozygous genotypes. The simulated values of the following pharmacokinetic parameters are outside the experimental range: t_{1/2, E-3174} for the CYP2C9*1/CYP2C9*2 group, C_{max, losartan}, AUC_losartan and t_{1/2, E-3174} for CYP2C9*1/CYP2C9*3, C_{max, E-3174} and t_{1/2, E-3174} for CYP2C9*2/CYP2C9*3. At the same time, the simulated AUC_E-3174 value for the CYP2C9*2/CYP2C9*3 genotype is at the upper bound of the mean ± SD experimental range.

To assess the antihypertensive effect of losartan therapy in individuals with different CYP2C9 genotypes, we used a previously developed cardiorenal model (Kutumova et al., 2022; Kutumova et al., 2021), which, in particular, reproduces the pharmacological action of losartan, and examined 100 virtual patients with arterial hypertension generated for it earlier (Kutumova et al., 2022). The distribution of their physiological characteristics is presented in Supplementary Figure S4. To test how virtual patients with different allelic variants of CYP2C9 would respond to losartan treatment, we estimated the values of the k _block parameter, representing the AT1-receptor blocking activity of the drug, for different CYP2C9 genotypes as described in Materials and Methods (Table 4). Each of the 100 virtual patients was simulated with each of the six values of k _block , i.e., the difference between the groups for each genotype consisted only in CYP2C9 activity.

After simulating the treatment of virtual patients with the found values of k _block , we compared our results with clinical study by Sinitsina et al. (2023), where authors had investigated the antihypertensive effect of losartan monotherapy in hypertensive patients with different CYP2C9 genotypes. In this study, patients were divided into two groups: 55 patients with the CYP2C9*1/CYP2C9*1 genotype and 26 patients carrying CYP2C9*2 or CYP2C9*3 variant alleles (13 patients with the CYP2C9*1/CYP2C9*2 genotype, 9 patients with the CYP2C9*1/CYP2C9*3 genotype, 2 patients with the CYP2C9*2/CYP2C9*2 genotype, and 2 patients with the CYP2C9*2/CYP2C9*3 genotype). Patients with the CYP2C9*3/CYP2C9*3 genotypes were not included in the study due to the rarity of this genotype.

Similar to the clinical study, we divided our virtual patients into the same two groups. For the CYP2C9*1/CYP2C9*1 group, we used all 100 patients with the corresponding k _block value (0.886), while for the second group, we used the same patients but with k _block values of other genotypes. Сomparison of the simulation and experimental results is presented in Table 7 and Figure 6. Unfortunately, we were unable to statistically evaluate the difference between our simulations and the clinical trial results because the experimental study (Sinitsina et al., 2023) only presented data as median and interquartile range. Overall, the agreement between modeled and experimental blood pressure responses to losartan treatment was better in patients with CYP2C9*2 or CYP2C9*3 alleles than in CYP2C9*1/CYP2C9*1 patients.

In addition, we compared the simulated systolic and diastolic blood pressure responses to losartan therapy between patients with the CYP2C9*1/CYP2C9*1 genotype and CYP2C9*2 or CYP2C9*3 carriers (Figure 7). As a result, we can conclude that virtual patients with the CYP2C9*1/CYP2C9*1 genotype had a greater reduction in both systolic and diastolic blood pressure after treatment with losartan, than patients with CYP2C9*2 or CYP2C9*3 alleles (p-value < 0.0001).

We then examined the antihypertensive effect of losartan across all CYP2C9 genotypes (Figure 8; Table 8). As can be seen, patients carrying the CYP2C9*2 allele had a similar reduction in blood pressure as patients with the CYP2C9*1/CYP2C9*1 genotype. At the same time, patients with the CYP2C9*3 allele had a significantly lower blood pressure response to losartan than patients with the wild-type genotype, while patients with CYP2C9*3/CYP2C9*3 did not experience any effect from losartan treatment.

The main limitation of this study is related to the small number of patients with different CYP2C9 genotypes whose clinical data were used for model validation and verification. This is because the frequencies of the variant alleles CYP2C9*2 and CYP2C9*3 are also very low.

Accumulating the results on modeling the influence of various genetic factors on the regulation of arterial pressure may further enable studying their cross-influence on the development and treatment of arterial hypertension. In particular, other genes that could influence the absorption, distribution, metabolism, and excretion of losartan can be considered, which will allow for more accurate prediction of the pharmacokinetics and pharmacodynamics of losartan.