A novel cardiovascular systems model to quantify drugs effects on the inter-relationship between contractility and other hemodynamic variables.

The use of systems-based pharmacological modeling approaches to characterize mode-of-action and concentration-effect relationships for drugs on specific hemodynamic variables has been demonstrated. Here, we (i) expand a previously developed hemodynamic system model through integration of cardiac output (CO) with contractility (CTR) using pressure-volume loop theory, and (ii) evaluate the contribution of CO data for identification of system-specific parameters, using atenolol as proof-of-concept drug. Previously collected experimental data was used to develop the systems model, and included measurements for heart rate (HR), CO, mean arterial pressure (MAP), and CTR after administration of atenolol (0.3-30 mg/kg) from three in vivo telemetry studies in conscious Beagle dogs. The developed cardiovascular (CVS)-contractility systems model adequately described the effect of atenolol on HR, CO, dP/dtmax, and MAP dynamics and allowed identification of both system- and drug-specific parameters with good precision. Model parameters were structurally identifiable, and the true mode of action can be identified properly. Omission of CO data did not lead to a significant change in parameter estimates compared to a model that included CO data. The newly developed CVS-contractility systems model characterizes short-term drug effects on CTR, CO, and other hemodynamic variables in an integrated and quantitative manner. When the baseline value of total peripheral resistance is predefined, CO data was not required to identify drug- and system-specific parameters. Confirmation of the consistency of system-specific parameters via inclusion of data for additional drugs and species is warranted. Ultimately, the developed model has the potential to be of relevance to support translational CVS safety studies.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Hemodynamic systems models in rat and dog have been established to characterize drug effects on the inter‐relationship between key hemodynamic biomarkers, such as blood pressure, heart rate, and cardiac output (CO). However, to date, no models exist that integrate contractility with hemodynamic variables in a systems manner.

WHAT QUESTION DID THIS STUDY ADDRESS?

Can the inter‐relationship between contractility and other hemodynamic variable be described in a systems manner? Can contractility readouts (dP/dtmax) replace CO readouts to identify drug‐ and system‐specific parameters?

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

A novel cardiovascular (CVS) systems model was developed based on pressure‐volume loop theory using atenolol as proof‐of‐concept drug. Both systems‐ and drug‐specific parameters of this CVS‐contractility (CTR) model could be identified without CO data. Hence, the absence of CO data does not limit the application of this model to data from other species.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

It is expected that the incorporation of pressure‐volume loop theory will provide the mechanistic basis to apply the developed CVS‐CTR model for inter‐species scaling. It has the potential to be integrated into a translational modeling platform to support CVS safety evaluations.

INTRODUCTION

Cardiovascular (CVS) safety issues are among the major causes of safety‐related attrition during drug development. , , In this context, most attention has been given to potentially fatal arrhythmias arising from QT‐prolongation, but non‐QT drug effects on hemodynamic end points, such as blood pressure or cardiac contractility, also represent serious CVS safety issues. , , , , , Preclinical in vivo telemetry studies are often conducted to identify both QT‐ and non‐QT effects, in which an ascending dose or placebo is administered to animals while measuring the time course of changes in heart rate (HR), blood pressure, cardiac output (CO), and/or contractility (CTR; e.g., dP/dtmax).

The interpretation of preclinical data from CVS safety studies, designed to determine the drug mode‐of‐action and concentration‐effect relationship for specific hemodynamic variables, can be challenging due to underlying complex homeostasis and feedback inter‐relationships. Quantitative CVS systems models to characterize the dynamics and pharmacokinetic/pharmacodynamic (PK/PD) relationships of hemodynamic variables after drug administration have previously been proposed to address this challenge. , A hemodynamic systems model was established for eight compounds with diverse modes‐of‐action in rats, here referred to as the “Snelder model.” , This model characterized drug effects on the inter‐relationship among HR, stroke volume (SV), total peripheral resistance (TPR), CO, and mean arterial pressure (MAP). Importantly, the Snelder model allows for separation of both drug‐ and system‐specific (e.g., drug‐independent) parameters. For identification of the system‐specific parameters, invasive and challenging CO measurements were required, which limits the integration of this model into a translational modeling platform, which would be an ultimate goal. Alternatively, measures for CTR, such as left ventricular (LV) dP/dtmax, could replace CO, as it can be measured more easily using telemetry. The dP/dtmax is determined by myocardial CTR and the loading conditions on the ventricle, and has been reported to be significantly correlated with SV and HR. Based on this, and expanding on the work from Snelder et al., an adapted model was recently published using data from telemetry studies in dogs, which replaced CO measurements by CTR (dP/dtmax) measurements. This model can be used to characterize novel compounds’ effects on the CVS system in dogs. However, in the absence of a mechanistic basis for the inter‐relationship among CO, SV, and CTR, the predictive value for other species could be limited.

LV pressure and volume have long been studied with pressure‐volume (PV) loop theory. A PV loop can be generated by plotting multiple measurements of LV pressure against LV volume in a complete cardiac cycle. In this theory, the changes in the end‐systolic pressure‐volume relationship reflect inotropic changes and the influence of vascular loading on systolic pressure and flow generation. This provides insight in the myocardial CTR and loading conditions of the heart for the hemodynamic modeling. As the PV loop theory captures the physiological principles of the cardiac cycle, it is uniquely suited as a basis for integration of CTR in the Snelder model in a system‐based manner.

Here, we describe the development of a novel CVS‐CTR systems model, which integrates dP/dtmax with CO, and other hemodynamic biomarkers (including HR and MAP) based on the principles from PV loop theory, using data from multiple dog telemetry studies for the selective β1‐blocker atenolol as proof‐of‐concept. We characterized the model through identifiability analyses with respect to drug‐ and system‐specific parameters, performed external validation studies, and investigated if dP/dtmax can replace CO for the identification of system‐specific parameters.

METHODS

Experimental data

We obtained data from three in vivo telemetry studies conducted at or on behalf of different pharmaceutical companies in conscious Beagle dogs, which were administered increasing doses of atenolol (0.3–30 mg/kg) or placebo with washout periods in between treatments (Table 1). Each study was conducted according to company‐specific standard procedures. In each study, several CVS hemodynamic variables were measured over time for a period of 6 to 24 h. In this analysis, we only used time course measurements for MAP, HR, CO, and LV dP/dtmax. Data from studies 1 and 2 were used for model development, whereas data from study 3 were used for external validation. The data from studies 1 and 3 were obtained after administration of control compounds as part of earlier drug development programs, and have not been previously published. The animals and experimental procedures are described for each respective study in the Supplementary Materials.

TABLE 1

Overview of studies used for model development

Study group	Site	Sex	n	Route of administration	Dose level	Dosing day	Sampling time range	Measurements
1	Servier	M	4	Oral	0 a , 3, 10, 30 mg/kg	1, 5, 7, 10	−1 to 6 h	HR, dP/dtmax, CO, MAP
2	AstraZeneca	M	4	Oral	0 b , 1, 3, 10 mg/kg	1, 5, 8, 10	−1 to 20 h	HR, dP/dtmax, MAP
3	GlaxoSmithKline	M	4	Oral	0 b , 0.3, 1, 3 mg/kg	1, 4, 7, 10	1 to 24 h	HR, dP/dtmax, MAP

Abbreviations: CO, cardiac output; HR, heart rate; MAP, mean arterial pressure.

Vehicle: Methylcellulose 0.5%.

Vehicle: water.

Pharmacokinetic model

Because PK was not measured in the experimental studies, we utilized a previously published PK model to predict drug plasma concentrations over time. The model was a three‐compartment model with first‐order absorption and elimination, , and the following PK parameters: Ka = 1.13 h−1, CL = 3.35 L/(kg*h), V2 = 4.05 L/kg, Q2 = 8.85 L/(kg*h), V3 = 3.74 L/kg, Q3 = 5.73 L/(kg*h), V4 = 11.9 L/kg, and F1 = 0.783. The predicted concentrations in the central compartment were directly linked to the drug effects.

Relationship between hemodynamic variables

We defined relationships between the hemodynamic variables used in the systems modeling framework as follows: MAP was described as a product of CO and TPR based on Ohm’s law, and CO was the product of HR and SV (Equation 1)

SV was defined further as:

To describe the inter‐relationship between SV and CTR, we considered the LV PV loop theory. There are two key linear relationships in PV loops (Figure 1). Line EA (Figure 1) represents the afterload relationship with the slope as arterial elastance (EA) (Equation 3) where EA is the product of HR and TPR (Equation 4).

FIGURE 1

Quantitative inter‐relationships in left ventricle (LV) pressure‐volume loop. (1) Red dashed line 1: the afterload relationship. The slope equals the arterial elastance (EA), which is the product of heart rate (HR) and total peripheral resistance (TPR). Its x‐axis intercept is end diastolic volume (EDV); (2) Blue dashed line 2: the end‐systolic pressure‐volume relationship (ESPVR) with the slope contractility (CTR) and x‐axis intercept of V0, which is the minimum volume needed to get pressure in left ventricle; PA: the pressure at the intersection point A of line 1 and line 2, when aortic valve closes in a single cardiac cycle; SV: stroke volume, the difference between EDV and end systolic volume (ESV)

Line ESPVR (Figure 1) represents the end‐systolic pressure‐volume relationship, which depends on CTR, and is defined as follows in relation to CTR:

Here, V0 is the x‐axis intercept of ESPVR, or the minimum volume required to get pressure in the LV. Equations 3 and 5 were solved to obtain the following relationship for ESV (Equation 6): dP/dtmax is used in this study as indicator of myocardial CTR, which depends on end diastolic volume (EDV; preload) and duration of isovolumetric contraction (TIC). , Therefore, we defined the variable CTR measurements (CTRM), which reflects for dP/dtmax measurements (Equation 7), as follows:

For a typical young Beagle dog (27 months) a TIC of 0.256 s was reported.

Model structure

The CVS system is regulated via various control mechanisms, such as the baroreflex system and the renin‐angiotensin system. To account for these negative feedback mechanisms, and to allow for the inclusion of drug effects into the model, the dynamics of HR, EDV, TPR, and CTR were described by a set of four differential equations, which were linked via negative feedback through MAP (Equation 8), similar to the Snelder model. ,

Here, FB is the feedback effect parameter through MAP to all the three primary parameters. The system was forced to be in steady‐state, or dynamic equilibrium, before drug administration, which means HR, EDV, TPR, and CTR do not change over time and equal their baseline (BSL) values. Therefore, the production rates of HR, EDV, TPR, and CTR can be expressed in terms of baseline and kout as follows:

Circadian rhythms in hemodynamic variables

Three cosine functions were defined as follows to model the circadian rhythm in HR, TPR, and CTR (CRHR, CRTPR, and CRCTR) influencing kinHR, kinTPR, and kinCTR, respectively:

Here, amp refers to amplitude and hor is the horizontal displacement. We tested periods of 8, 12, and 24 h, and identified the optimal period for each variable. The system was forced to be in oscillating steady‐state at the start of pharmacological intervention by shifting the observations 1 week (determined empirically) in time (i.e., the system was initialized at 0 h and the pharmacological interventions started at 168 h).

Modelling of the concentration‐effect relationship

Concentration‐effect relationships were evaluated using maximum effect (Emax) models.

Here, EFF represents the drug effect at concentration C, Emax is the maximum effect constant, and EC50 is the concentration constant at which half of the maximum effect is achieved. Based on literature values for the binding affinity (KD) of atenolol to the β1‐ and β2‐adrenoceptors from Baker et al., we fixed the EC50 for the effects on HR and CTR, which are mediated through the β1‐adrenoceptors, to 58.3 ng/ml and the EC50 for the effect on TPR, which is mediated through the β2‐adrenoceptor, to 272.5 ng/ml.

Incorporation of drug effects in the systems model

The complete model consisted of a combination of Equations (1), (2), (3), (4), (5), (6), (7), (8), (9), (10), (11) and the following differential equations:

Structural identifiability analysis

A structural identifiability analysis was performed to verify whether the model parameters could be estimated from experimental data. Structural identifiability of the model was evaluated using the MATLAB toolbox GenSSI 2.0, which is a MATLAB toolbox together with MATLAB (version R2020a). We implemented simultaneously two different drug effects, namely, EFFHR and EFFCTR, to influence the production rate kin of HR and CTR, respectively. At the same time, we applied three different circadian rhythms (i.e., CRHR, CRTPR, CRCTR), to influence the production rate of HR, TPR, and CTR, respectively. We verified the structural identifiability of the model for two different sets of observable outputs, namely HR, CTRM, CO, and MAP and HR, CTRM, and MAP.

Model estimation and development

The data from studies 1 and 2 were simultaneously analyzed using a nonlinear mixed‐effects modeling approach implemented in NONMEM (version 7.4.3, Icon Development Solutions, Ellicott, MD) with Perl speak to NONMEM toolkit (PsN, version 4.8.1, Uppsala University, Sweden). First order linear conditional estimation method with interaction was used for model estimation. Interindividual variability (IIV) on BSLHR, BSLTPR and BSLCTRM were evaluated using log‐normal distributions. Exponential, additional, proportional, and combined error models were evaluated for residual error variability. A decrease of 3.84 (corresponding to p < 0.05 in a chi‐squared distribution with degree of freedom = 1) in the objective function value (OFV) by adding an additional parameter between nested models was considered statistically significant. The baseline parameters for HR, CTRM, and TPR were estimated together with V0, while ensuring that these parameters remain positive during the estimation.

Relevance of the CO data

To investigate the importance of CO data in the estimation, the final model was fitted to a dataset without CO readouts. In the absence of remaining informative data, the circadian rhythm (CRTPR) and IIV on BSLTPR were omitted from the model and BSLTPR was fixed to 0.0743 mmHg*ml/min, which was the estimate in the final model. Fixing this parameter was required to avoid over‐parameterization.

Model evaluation

Visual predictive checks (VPC) with 1000 samples and goodness‐of‐fit plots (GOF) were conducted to evaluate the predictive performance of the final model. External validation consisted of comparing the model prediction (95% prediction interval using the final model) with hemodynamic observations from study 3.

External validation

The external validation was conducted based on the data from study 3. A VPC was performed to predict the effect of atenolol on the hemodynamic markers as detailed here above. All the system‐specific parameters were assumed the same as the final estimates from study 1.

Identification of the right site of action of compounds given the data

Stochastic simulations and re‐estimation (SSE) using PsN toolkit in conjunction with NONMEM were conducted to investigate if the model can be used to identify the mode of action of new drugs. We simulated three types of models with drug effects on HR, CTR, or TPR, defined as the original models in this analysis. Subsequently, we re‐estimated the parameters on the simulated datasets using models with drug effects on HR, CTR, or TPR and no drug effect, defined as the alternative models in this analysis. The simulation scenarios included three combinations of measurements at rich timepoints for (i) HR, dPdtmax, CO, and MAP; (ii) HR, dPdtmax, and MAP; and (iii) HR and MAP. The parameter values used in the simulation steps were fixed to the final estimates from previous model fitting. The difference in OFV (dOFV) between drug effect models and no drug effect models was calculated to compare the performance of the models with different sites of action.

Simulations to illustrate the properties of the CVS system

To illustrate the properties of the system for drugs with a primary effect on HR, CTR, or TPR, simulations were performed using the final model (Equations (1), (2), (3), (4), (5), (6), (7), (8), (9), (10), (11)). The values of system‐specific parameters in the simulation were fixed to the final estimates obtained from the final model, and hypotheses of both negative and positive effects (reflecting different mechanisms of action [MOA]) were tested.

RESULTS

The structural identifiability analysis showed that the proposed model is at least locally identifiable for both sets of observable outputs (HR, CTRM, CO, and MAP and HR, CTRM, and MAP), indicating the possibility of estimating parameter values from experimental data. The difference between the two sets of observable variables is that when CO is measured all baseline parameters (i.e., BSLCTRM, BSLHR, BSLTPR, and V0) are globally identifiable indicating that they can be uniquely determined, whereas when no CO data are available BSLTPR and V0 are only locally identifiable indicating that more than one solution is possible.

Model development and evaluation

Data from studies 1 and 2 were adequately described by the final CVS‐CTR model (Figure 2 and Equations (1), (2), (3), (4), (5), (6), (7), (8), (9), (10), (11)) as can be seen in the VPC (Figure 3) and GOF plots (Figures S1–S4). All parameter could be estimated precisely, with relative standard errors less than 50% (Table 2). The concentration‐effect relationships were best described by Emax models. The EC50 for HR and CTR was fixed to the reported KD of 58.3 ng/ml for binding of atenolol to the β1 receptor. Emax for the effect on HR and CTR was estimated to be 0.415 and 0.422, respectively. An Emax model was added for a positive effect on TPR with the EC50 fixed to 272.5 ng/ml. However, as the final estimate of Emax for TPR approached zero, and the OFV did not decrease significantly (3.84, p < 0.05), this additional drug effect was rejected. Overall, the final model was structured with inhibiting effects on both HR and CTR.

FIGURE 2

CVS‐contractility model structure to characterize drug effects on the cardiovascular systems. HR, heart rate; CTR, contractility; CTRM, contractility measurements variable for dP/dtmax; TPR, total peripheral resistance; MAP, mean arterial pressure; CO, cardiac output; SV, stroke volume; ESV, end‐systolic volume; EDV, end‐diastolic volume; EA, arterial elastance. HR, TPR, and CTR are regulated by negative feedback through MAP, where FB represents the magnitude of feedback on HR, TPR, and CTR. Effects on HR, CTR, EDV, and TPR are described by four linked turnover equations, in which , , , and represent the zero‐order production rate constants, and , , , and represent the first‐order degradation rate constants of HR, EDV, TPR, and CTR, respectively. CTRM is the variable for dP/dtmax, which is influenced by CTR and EDV. The inter‐relationships between ESV, EDV, EA, and CTR are illustrated in Figure 1. Potential drug effect cites are indicated in orange. Hemodynamic variables which can be observed are indicated in light blue

FIGURE 3

Visual predictive check for the final model. (a) Study 1 data following administration of placebo, 3 mg/kg, 10 mg/kg, and 30 mg/kg atenolol. (b) Study 2 data following administration of placebo, 1 mg/kg, 3 mg/kg, and 10 mg/kg atenolol. The grey points and grey lines represent the observed data for biomarkers of HR, dP/dtmax, CO, and MAP. The black lines represent the median of the observed data for each biomarker. The shaded blue areas represent the 95% CI of the median of the predictions. CI, confidence interval; CO, cardiac output; HR, heart rate; MAP, mean arterial pressure

TABLE 2

The system‐ and drug‐specific parameter values from the CVS‐contractility model fitted to a dataset including CO data (final model) and a dataset from which the CO data were omitted (model without CO data)

Parameter		Final model		Model without CO data
		Value (RSE)	95% CI	Value (RSE)
System‐specific parameters
BSLHR_1 (beats*min−1)	Baseline of HR for study 1	79.4 (5.78%)	70.4–88.3	78.4 (4.69%)
V0_1 (ml)	V0 for study 1	9.92 (9.81%)	8.02–11.8	10.1 (11.0%)
BSLCTRM_1 (mmHg/s)	Baseline of dP/dtmax for study 1	3777 (7.59%)	3215–4338	3652 (9.04%)
BSLHR_2 (beats*min−1)	Baseline of HR for study 2	77.0 (2.63%)	73.0–81.0	76.4 (2.73%)
V0_2 (ml)	V0 for study 2	9.15 (9.62%)	7.42–10.9	8.55 (4.76%)
BSLCTRM_2 (mmHg/s)	Baseline of dP/dtmax for study 2	2422 (7.50%)	2067–2779	2480 (7.70%)
BSLTPR (mmHg*min/ml)	Baseline of TPR	0.0743 (4.55%)	0.0677–0.0810	0.0743 FIX
BSLEDV (ml)	Baseline of EDV	31.13	FIX	31.13 FIX
Kout a (h−1)	Degradation rate	0.830 (21.8%)	0.475–1.18	1.07 (22.9%)
FB (mmHg−1)	Feedback effect	0.00558 (12.5%)	0.00422–0.00694	0.00589 (27.9%)
Circadian Rhythm
Amp_1 b	Amplitude for study 1	0.0931 (41.3%)	0.0178–0.168	0.0775 (48.1%)
HorHR_1	Horizontal displacement of HR	7.86 (28.7%)	3.43–12.3	6.33 (17.8%)
HorCTR_1	Horizontal displacement of CTR	9.82 (18.9%)	6.19–13.5	8.57 (16.6%)
Amp_2	Amplitude for study 2	0.168 (13.6%)	0.123–0.213	0.172 (22.6%)
HorHR_2	Horizontal displacement of HR	19.4 (1.60%)	18.8–20.1	19.5 (1.99%)
HorCTR_2	Horizontal displacement of CTR	21.8 (1.84%)	21.0–22.6	21.8 (2.84%)
HorTPR	Horizontal displacement of TPR	6.33 (1.96%)	6.09–6.58	0 FIX
Drug‐specific parameters
EC50 HR and EC50 CTR (ng/ml)	EC50 on HR and CTR	58.3	FIX	FIX
EmaxHR	Emax on HR	0.415 (11.6%)	0.321–0.510	0.402 (18.0%)
EmaxCTR	Emax on CTR	0.422 (9.56%)	0.343–0.501	0.421 (14.8%)
Interindividual variability
BSLHR (CV%)		6.08 (24.2%)	4.41–7.39	3.89 (45.0%)
BSLCTRM (CV%)		16.58 (34.2%)	9.47–21.52	16.99 (33.5%)
BSLTPR (CV%)		6.48 (29.6%)	4.2–8.15	0 FIX
BSLCTRM X BSLTPR (CV%)		8.88 (30.9%)	5.57–11.26	0 FIX
Residual variability
Res. ErrorHR (CV%)		11.53 (12.7%)	9.98–12.9	11.37 (11.9%)
Res. ErrordPdtmax (CV%)		10.43 (13.6%)	8.93–11.74	10.26 (13.6%)

Abbreviations: CI, confidence interval; CO, cardiac output; CVS, cardiovascular; EC50, half‐maximal effective concentration; EDV, end diastolic volume; Emax, maximum effect; HR, heart rate; RSE, relative standard error; TPR, total peripheral resistance.

KoutHR = KoutCTR = KoutTPR = KoutEDV = Kout.

AmpHR = AmpCTR = AmpTPR = Amp, Amp_1 and Amp_2 are the amplitudes of circadian rhythm for studies 1 and 2 respectively, and for TPR the Amp and Hor were assumed to be the same as in study 1.

The rate constants KoutHR, KoutCTR, KoutTPR, and KoutEDV were initially estimated separately and did not differ significantly. As the OFV did not increase significantly (p < 0.05) after estimating the same Kout for all four biomarkers, these four parameters were lumped into one parameter (Kout). Because CO was not measured in study 2, the typical values of BSLTPR were assumed valid for both studies. The other baseline parameters differed between the two studies, with a 3.12% higher BSLHR and 8.42% higher V0 for study 1 as compared to study 2, whereas BSLCTRM for study 1 was 55.94% higher than for study 2. BSLSV, BSLCO, and BSLMAP were derived from these parameters, resulting in higher BSLSV, BSLCO, and BSLMAP for study 1 as compared with study 2. The amplitude parameters ampHR, ampCTR, and ampTPR, could not be distinguished and were estimated to be 0.0931 for study 1 and 0.168 for study 2.

The IIV was estimated on the baseline parameters BSLHR, V0, BSLCTRM, and BSLTPR, and covariance was estimated between ETA_BSLHR and ETA_V0, and between ETA_BSLCTRM and ETA_BSLTPR. The residual errors on HR and dP/dtmax were best described by exponential residual error models with lower RSEs compared to additive, proportional, and combined residual error models. The residual error of TPR was very small and, therefore, it was fixed to 0.

Evaluation of the relevance of the CO data

Fitting the CVS‐CTR model to a dataset from which the CO data were omitted did not result in a significant change in parameter estimates; almost all estimates were within the 95% confidence interval (CI) of the parameters of the final model fitted to the complete dataset including CO data (Table 2). The only parameter outside the 95% CI, IIV on BSLHR, was also close to the lower boundary of 95% CI and would not significantly influence the model prediction.

The external validation for study 3 (Figure S5) showed that in general the validatory data were reasonably well‐predicted. However, the model slightly underpredicted the effect on dP/dtmax in the 3 mg/kg dose group, possibly reflecting interstudy variability.

In the SSE analysis, the models with the “correct” MOA showed the highest reduction in OFV for all three simulation scenarios (Figure 4). With decreasing numbers of observations, the difference of “correct” and “incorrect” models in OFV gradually shrank. Especially, when re‐estimating the parameters on the HR/MAP data from the model with the original effect on CTR, the dOFV for “incorrect” and “correct” models was quite close.

FIGURE 4

Identification of the site of action of compounds with effects on HR, CTR, or TPR. The y‐axis shows the delta objective function value of the re‐estimated models compared to a model with no drug effect. The red bars are incorrect alternative models, which do not match the original model, and the green bars are correct alternative models, which match the original model. The SD bars represent the standard error around the mean value. CTR, contractility; HR, heart rate; MAP, mean arterial pressure; TPR, total peripheral resistance

In the simulations of drug effects with different MOAs, there were noticeable differences between the signature profiles of HR, dP/dtmax, CO, TPR, and MAP. As shown in Figure 5, inhibition of HR, CTR, or TPR always lead to a drop in MAP indicating that the magnitude of negative feedback from MAP is lower than the primary effect. Similarly, the positive effects on HR, CTR, or TPR resulted in an increase in MAP, as shown in Figure S6. The site of action of a compound with unknown MOA can be derived from the combination of the direction of the effects on MAP, HR, and CTR.

FIGURE 5

Signature profiles following administration of 30 mg of a hypothetical compound with negative effect on HR (panel a), CTR (panel b), and TPR (panel c), respectively. The red lines represent a decrease the biomarkers after pharmacological intervention. The green lines represent an increase in the biomarkers after pharmacological intervention. CO, cardiac output; CTR, contractility; HR, heart rate; MAP, mean arterial pressure; TPR, total peripheral resistance

DISCUSSION

We describe the development of a novel CVS‐CTR systems model to characterize drug effects on the inter‐relationship among CTR and CO, MAP, HR, TPR, SV, ESV, and EDV. The proposed model was based on a previously developed CVS model by Snelder et al., , which captures the inter‐relationship between CO, MAP, HR, TPR, and SV. To date, no mechanistic models exist that integrate CTR with the other mentioned hemodynamic variables, except for a model that was recently published by Venkatasubramanian et al. This model describes SV as a function of CTR and HR, to account for the effect of lowering SV due to shorter filling time with increased HR during diastole, and was successfully applied to characterize drug effects on the hemodynamic system in dogs. Although the relationship between SV, CTR, and HR has a mechanistic basis, the proposed model is not truly mechanistic, as it does not account for the effect of afterload/TPR on the ability of the ventricle to eject blood, altering ESV and SV. PV‐loop theory does cover these aspects. Hence, we based our model on this theory. Better capturing the physiology of the CVS system via PV‐loop theory allows for improved separation of system‐specific and drug‐specific parameters, which is crucial to understand variation between species. Therefore, the developed CVS‐CTR model is the first systems model that integrates CTR with other key hemodynamic variables. By fitting the CVS‐CTR systems model to a dataset that did not include CO data, we showed that CTR data can replace CO data. The absence of CO data does not limit the application of this model to data from other species, which facilitates the development of a translational modeling platform.

The developed CVS‐CTR model adequately described the effect of atenolol on the inter‐relationship between CTR and the other hemodynamic variables. Atenolol was selected as a proof‐of‐concept compound because it is widely used, and its hemodynamic effects are well‐investigated. To identify the system‐specific parameters in the best possible way, we fixed the EC50 for the effects of atenolol on HR and CTR to literature (KD) values. Although the EC50 is a lumped parameter for binding to the receptor and signal transduction, we assumed that the EC50 could be fixed to the reported KD. To evaluate if these KDs were covered by the exposures in the dog studies we simulated the atenolol PK profiles following single oral doses of 0.3, 1, 3, 10, or 30 mg/kg atenolol (Figure S7). Simulated concentrations after doses of 0.3 and 1 mg/kg did not reach the KD of β1‐adrenoceptors, which is in line with the weak drug effect observed in these low‐dose groups. Moreover, at high doses of 10 and 30 mg/kg, only peak concentrations were above the KD for binding to β2‐adrenoceptors, which is in line with the finding that adding a positive drug effect on TPR did not significantly improve the description of the data. Overall, our assumed EC50 seems justified, as the CVS effects of atenolol could be adequately described by a combination of primary inhibitory effects on both HR and CTR, which corresponds well with its high affinity for β1‐adrenoceptors. By comparing the value of amplitude and Emax parameters, it was found that the magnitude of circadian rhythm is approximately three to four times lower than the maximum drug effect. In addition, in the current analysis we have fixed KoutCTR, KoutTPR, and KoutEDV to the same value as KoutHR, as these parameters could not be estimated independently. Identical Kout values could be biologically plausible because all variables studied are influenced by feedback through the baroreflex system. To confirm this, our model should be fitted to data for additional drugs with diverse mode of action in the future.

In our experimental data, a drop was observed in both HR and CTR after administration of atenolol. Therefore, we assumed that the primary drug effect is on both HR and CTR, which is in line with the MOA of β1‐blockers. This assumption was challenged by forward inclusion and backward exclusion during the model estimation. An increase in HR leads to an increase in the force of contraction generated by the myocardial cells through the Bowditch effect, which indicates that CTR can be regulated by HR. This phenomenon is associated with calcium ion handling and mishandling in cardiac cells. However, this pathway was not included in our current model structure as it could not be distinguished from the effect atenolol on both HR and CTR. In future research, we will consider to introduce this regulation on CTR by HR in the model, if a similar behavior pattern of HR and CTR is observed for other compounds with different MOAs as well.

The developed model incorporates EDV as a key variable. As the available data did not contain any information to estimate baseline EDV, this parameter had to be fixed to a reported literature value to avoid over‐parameterization. Importantly, although EDV does not currently impact model predictions, its inclusion was considered necessary to allow future application, for example, compounds with a primary drug effect on EDV. Dynamics of EDV may be further evaluated in future work if data for drugs with a known effect on SV would be available, and the necessity to estimate a different Kout for EDV could then be investigated.

In the MOA‐identification analyses, the site of action of new compounds could be identified correctly with three combinations of observations: (1) HR, CO, dP/dtmax, and MAP; (2) HR, dP/dtmax, and MAP; and (3) HR and MAP. This implies that CO readouts are not needed for identification of the site of action of compounds in the SSE analyses. In the simulations to illustrate the system properties after simulating drug effects on HR, CTR, and TPR, there are noticeable differences in the behavior of HR, dP/dtmax, TPR, CO, and MAP (Figure 5 and Figure S5). This indicates that the site of action of new compounds can be distinguished by comparing the direction of effect on HR, dP/dtmax, and MAP. Because the negative feedback effect is always lower than the primary effect, the primary site of action is the biomarker with the same direction of effect as the effect on MAP. Thus, we anticipate that the model can be used to identify the site of new compounds with a single site of action, for which the mechanism underlying CVS effect is unclear. For more complex MOAs with more than one primary site of action, different assumptions will need to be evaluated.

CONCLUSION

We successfully developed a novel CVS‐CTR model to characterize drug effects on the inter‐relationship between CTR and other key hemodynamic variables based on pressure‐volume loop theory. Through fitting our model to data from telemetry studies of experimental compounds, our model can support the identification and quantification of the mode of action on key hemodynamic end points. We expect that the incorporation of PV‐loop theory will provide mechanistic basis that could be of relevance when applying this model for interspecies scaling. Key hurdles in this extrapolation may, for instance, be differences in homeostatic feedback mechanisms and set points. Future studies may focus on applying the developed model to additional compounds in different species to confirm consistency of system‐specific parameters and scalability toward other species, and to potentially study how disease state or exercise will influence hemodynamic responses.

CONFLICT OF INTEREST

The authors declared no competing interests for this work.

AUTHOR CONTRIBUTIONS

Y.F., H.T., M.M.S., J.G.C.H., and N.S. designed and performed the research and analyzed the data. Y.F., H.T., M.M.S., E.R., T.A.C., S.B., D.L., P.H.G., J.G.C.H., and N.S. wrote the manuscript.

Supplementary Material

Click here for additional data file.