Abstract
Cardiovascular adverse effects in drug development are a major source of compound attrition. Characterization of blood pressure (BP), heart rate (HR), stroke volume (SV), and QT-interval prolongation are therefore necessary in early discovery. It is, however, common practice to analyze these effects independently of each other. High-resolution time courses are collected via telemetric techniques, but only low-resolution data are analyzed and reported. This ignores codependencies among responses (HR, BP, SV, and QT-interval) and separation of system (turnover properties) and drug-specific properties (potencies, efficacies). An analysis of drug exposure–time and high-resolution response-time data of HR and mean arterial blood pressure was performed after acute oral dosing of ivabradine, sildenafil, dofetilide, and pimobendan in Han-Wistar rats. All data were modeled jointly, including different compounds and exposure and response time courses, using a nonlinear mixed-effects approach. Estimated fractional turnover rates [h-1, relative standard error (%RSE) within parentheses] were 9.45 (15), 30.7 (7.8), 3.8 (13), and 0.115 (1.7) for QT, HR, total peripheral resistance, and SV, respectively. Potencies (nM, %RSE within parentheses) were IC50 = 475 (11), IC50 = 4.01 (5.4), EC50 = 50.6 (93), and IC50 = 47.8 (16), and efficacies (%RSE within parentheses) were Imax = 0.944 (1.7), Imax = 1.00 (1.3), Emax = 0.195 (9.9), and Imax = 0.745 (4.6) for ivabradine, sildenafil, dofetilide, and pimobendan. Hill parameters were estimated with good precision and below unity, indicating a shallow concentration-response relationship. An equilibrium concentration–biomarker response relationship was predicted and displayed graphically. This analysis demonstrates the utility of a model-based approach integrating data from different studies and compounds for refined preclinical safety margin assessment.
SIGNIFICANCE STATEMENT A model-based approach was proposed utilizing biomarker data on heart rate, blood pressure, and QT-interval. A pharmacodynamic model was developed to improve assessment of high-resolution telemetric cardiovascular safety data driven by different drugs (ivabradine, sildenafil, dofetilide, and pimobondan), wherein system- (turnover rates) and drug-specific parameters (e.g., potencies and efficacies) were sought. The model-predicted equilibrium concentration–biomarker response relationships and was used for safety assessment (predictions of 20% effective concentration, for example) of heart rate, blood pressure, and QT-interval.
Introduction
In drug development, cardiovascular effects are one of the major sources of compound attrition. Identifying these effects at an early stage is therefore a focus of significant effort (Piccini et al., 2009; Gintant et al., 2016). Preclinical data on cardiovascular variables, such as blood pressure (BP), heart rate (HR), and QT-interval prolongation, have long been routinely collected in animals and humans (Piccini et al., 2009). Although these variables are often obtained simultaneously, they are commonly analyzed separately. However, since significant codependencies between safety variables are typically shown, simultaneous modeling of all data results in more informative insights about the test substance (e.g., potency EC50) and systems parameters [e.g., fractional turnover rates (kout)]. Codependencies act both directly, as in QT adaptation to HR change, and indirectly (e.g., via baroreflex feedback) (Cleophas, 1998; Holzgrefe et al., 2014). In addition, useful information about primary drug effects is neglected when codependencies are ignored. A consequence of this is inaccurate results. It is also common practice to analyze the outcomes of experiments independently from each other in spite of repeatedly using the same test population, which partly ignores valuable baseline behavior. Furthermore, high-resolution time courses of biomarker data are collected using telemetry, but only low-resolution data are analyzed and reported. To remedy some of these shortcomings, pharmacokinetic (PK)/pharmacodynamic modeling has proven to be a valuable tool combining quantities, such as HR, BP, and QT, using a nonlinear mixed-effects (NLME) approach.
Semimechanistic models of cardiovascular physiology have been used for over 60 years, beginning with pivotal work by Noble (1962) and Guyton et al. (1972). Their work has been developed in a number of ways in response to improved understanding, new measurement techniques, and computational power (Kappel and Peer, 1993; Winslow et al., 1999; Ursino and Magosso, 2003; ten Tusscher, 2004; ). Models that accommodate drug effects have been introduced more recently (Zemzemi et al., 2013; Snelder et al., 2014; Colatsky et al., 2016). The present analysis extends the pivotal Guyton et al. (1972) framework used within many discovery projects and emphasizes the applicability of presently used models with high-resolution data. The ambition is not to use the collected data for model discrimination purposes but rather to demonstrate how a semimechanistic platform well known among physiologists and (safety) pharmacologists can be used with their data.
The model-based approach for improved assessment of cardiovascular safety that we propose here is shown schematically in Fig. 1. In step 1, high-resolution telemetric data are collected from Han-Wistar rats, and exploratory and regression analyses are done to visualize and quantify baseline properties of biomarker data, such as HR, BP, and QT-interval (Kramer and Kinter, 2003). In step 2, low-resolution biomarker data obtained with different drugs (ivabradine, sildenafil, dofetilide, and pimobendan) that have different mechanisms of action are analyzed using all biomarker data and different tool compounds simultaneously, thereby yielding separate estimates of drug and system properties. In step 3, equilibrium predictions of concentration–biomarker response relationships are calculated, visualized graphically, and quantified (e.g., in terms of model-derived EC20 for drug effects on HR, BP, and QT-interval). This step serves as the basis for future safety assessment.
Schematic overview of study analyses. Step 1: High-resolution time courses of biomarker response are obtained telemetrically from animals. These data are used to assess baseline biomarker response. Step 2: Low-resolution data combined with baseline model serve as input to a meta-analysis of biomarker-response time courses, compounds, and experiments, generating final estimates of systems (e.g., turnover rates) and drug (e.g., potencies) parameters. Step 3: Steady-state concentration–biomarker response relationships are established for heart rate, blood pressure, and QT response.
Materials and Methods
Animals
Male rats (Han-Wistar, Charles River) weighing 450–600 g, 3–19 months old, were used. Rats were chronically instrumented for the telemetric collection of arterial pressure, ECG, and body temperature as described previously (Schierok et al., 2000). Groups of up to three rats were housed together in cages (type IV—1815 cm2) in a room with a 12-hour light/dark cycle and controlled temperature and humidity, which was identical for all experiments. The rats had access to normal rodent chow and water ad libitum and were not fasted prior to experiments.
Telemetric System
The telemetric system consisted of elements of a Data Science International (DSI) system (DSI, Minnesota, USA) for cardiovascular measurements (BP, HR, temperature, and ECG). A calibration amplifier was added that was capable of converting the signals from the DSI system for transfer to a computer-based acquisition and analysis system (HEM v.3.3 or v.4.3, Notocord, Paris, France). Two of these systems were run in parallel to allow simultaneous measurement of eight animals.
Miniature transmitters (DSI TA PA-C40) for the measurement of BP and HR, temperature, and ECG were used. Before implantation, the transmitter was calibrated using an assigned receiver (RLA 1020 DSI). Afterward, the transmitters were sterilized with ethylene oxide, and after 24–48 hours of venting, they were kept under aseptic conditions until implantation.
Implant Surgery
Rats were anesthetized in a box ventilated with 4% isoflurane. After the onset of the anesthesia, they were transferred to the head chamber of the inhalation apparatus for anesthesia in small animals, Narcoquip (Völcker, Germany), and anesthesia was maintained with 2% isofluorane. Under aseptic conditions, a midline abdominal incision (2 and 3 cm) along the linea alba was made, and the intestines were retracted. The lower abdominal aorta was isolated and temporarily occluded with a ligature. A small hole was punctured into the aorta near the iliac bifurcation using a 20-gauge hypodermic needle, and the catheter end of the transmitter was inserted. After the puncture site was dried thoroughly, the catheter was sealed in place using a tissue adhesive (Histoacryl, Braun, Germany). The intestines were gently put back in place, and the transmitter body was attached to the abdominal muscle during closing of the abdominal cavity and the skin incision. To reduce any pain or infection risk, the animals were treated with dypyrone 50 mg/kg i.p. (Novalgin, Hoechst, Germany) prior to the surgery, and a prophylactic dose of penicillin 0.05 ml/100 g (Tardomyocel, Bayer, Germany) was administered. The surgical preparation was followed by a 4-week recovery period before experiments.
Study Description
Data from four separate drug provocation studies in rats were used: ivabradine (Sigma-Aldrich, US), sildenafil (Neuraxpharm, Germany), dofetilide (Sigma-Aldrich), and pimobendan (Sigma-Aldrich). Each compound was used in a pharmacokinetic and a pharmacodynamic study. We are consistently reporting and utilizing unbound plasma concentrations throughout this work, with unbound fractions listed in Table 1.
Compound, mode of action, doses, study period, molecular weight, and unbound fraction in rats
Unbound fractions from in-house data and FDA approval documentation (www.pharmapendium.com).
Pharmacokinetic Assessment
The plasma PK of all four compounds was investigated in rats after intravenous and oral administration. A single dose of test compound was either administered orally (volume administered: 2 ml/100 g b.wt. suspended in Natrosol 0.5% + 0.015% Tween80, n = 3 subjects per study) or intravenously (volume administered: 0.5 ml/100 g b.wt.; dissolved in a 20% HP-β-cyclodextrin solution with pH adjustment to pH 6, n = 2 subjects per study). For ivabradine, sildenafil, dofetilide, and pimobendan, doses were 1, 10, 1, and 3 mg/kg (oral) and 0.5, 5, 0.5, and 0.33 mg/kg (intravenous), respectively. Plasma drug concentrations were analyzed in samples drawn at 0.08, 0.25, 0.5, 1, 2, 4, 8, and 24 hours. Lower limits of quantification were 1 nM for ivabradine, dofetilide, and pimobendan and 2.5 nM for sildenafil.
Blood samples were collected in Microvette tubes 0.50 ml K3EDTA (Sarstedt). The tubes were ice-chilled prior to sampling. After collection, the blood was gently mixed by inverting the tube several times and stored upright on ice until centrifugation. Blood samples were centrifuged for 5 minutes at 8300g at 4°C. After centrifugation, plasma was aliquoted and stored below −20°C pending chemical analyses. Drug plasma concentrations of all four compounds were analyzed using liquid chromatography/tandem mass spectrometry with an analytical range of 1–2000 nM.
Pharmacodynamic Assessment
Four separate groups of rats (n = 8 per group) received vehicle control (Natrosol 0.5% + 0.015% Tween80) and one of the test compounds—ivabradine, dofetilide, pimobendan, or sildenafil—at three dose levels (volume administered: 1 ml/100 g b.wt. suspended in Natrosol 0.5% + 0.015% Tween80) orally (Table 1). Each study comprised eight individual animals except for the sildenafil study, which comprised 31 individual animals.
Eight implanted animals were transferred to eight separate cages placed on the previously assigned DSI receivers on the day of each experiment. The transmitters were activated magnetically, and the signal was confirmed with the help of a commercially available radio receiver. Systolic blood pressure (sBP), diastolic blood pressure (dBP), and ECG were continuously measured. From these data, three cardiovascular variables were derived: HR, mean arterial BP (MAP) (derived as one-third sBP + two-thirds dBP), and the QT-interval. All measurements began 1.5 hours before dosing and continued for variable lengths of recording (Table 1).
Data Acquisition
The Acquisition Software HEM version 3.3 or 4.3 (Notocord, Paris, France) was used to record the experimental data continuously (sampled at 1000 Hz) in real time and store it on the local hard disk. HEM is used to acquire hemodynamic arterial BP signals as well as ECG signals. From the arterial BP signal, systolic, diastolic, and mean arterial BP was calculated as well as HR. These signals are denoted high-resolution biomarkers. Body temperature was recorded directly. At the end of each experiment, the raw data were saved on a system server. The data were summarized by calculating the median value of sequential events over 10 minutes for each biomarker, giving low-resolution biomarkers.
Pharmacodynamic Modeling
Drug Exposure Model.
A schematic plot of the absorption and disposition models of ivabradine, sildenafil, dofetilide, and pimobendan is shown in Fig. 2.
Schematic diagram of the extended model. Left plot: absorption and disposition model used for plasma exposure-time data of ivabradine, sildenafil, and dofetilide. Right plot: absorption and disposition model used for plasma exposure-time data of pimobendan. For pimobendan, a hypothetical metabolite compartment (eq. 1d) was added to capture active metabolite(s) from parent compound. The rate constant kmet was estimated from biomarker response–time data and is not part of the plasma exposure analysis (Kitzen et al., 1988).
Drugs (also called test compounds) were administered intravenously or orally, and a two-compartment disposition model with first-order input/output was fitted to plasma concentration data:
(1a)
(1b)
(1c)
(1d)Ag, Cp, and Ct denote amount in the gut compartment, concentration in plasma, and concentration in tissue, respectively. The model parameters correspond to absorption rate constant (Ka), plasma clearance (CL), oral bioavailability (F), plasma volume of distribution (Vp), intercompartmental distribution (CLd), and tissue volume (Vt).
Note that a hypothetical metabolite compartment (eq. 1d) was added to capture active metabolite(s) from the pimobendan parent compound. The rate constant kmet was estimated from response-time data and is not a part of the plasma exposure analysis (Fasanmade and Jusko, 1995).
Pharmacodynamic Model
A schematic plot of the pharmacodynamic actions of ivabradine, sildenafil, dofetilide, and pimobendan is shown in Fig. 3. The model consists of four biomarker responses: HR, TPR, stroke volume (SV), and QT. HR, TPR, and SV are interconnected and drive MAP, which attenuates the turnover rates of these states, creating a feedback loop. HR also affects the turnover rate of the QT-interval, describing the length adaptation of QT to HR (Fig. 3).
Schematic diagram of the extended model. HR, SV, TPR, and QT-interval (QT) are described with their respective zero-order turnover rates (e.g., kin,HR) and first-order fractional turnover rates (e.g., kout,HR). The experimentally measured biomarker responses, such as heart rate, mean arterial pressure, and QT-interval are shown as dark gray boxes. There are no experimental data on total peripheral resistance and stroke volume (light gray boxes). Mean arterial pressure is also a derived function of HR, SV, and TPR. MAP acts on HR, SV, and TPR via negative feedback (dashed lines). The QT-interval is inhibited by HR (dashed line). Drug action on HR, SV, TPR and/or QT-interval is either via stimulatory action (e.g., S(P), in which pimobendan stimulates the turnover rate of SV) or via inhibitory action (e.g., I(I), in which ivabradine inhibits turnover rate of HR). Sildenafil and pimobendan both exhibit inhibitory action turnover rate of TPR (I(S, P)). Dofetilide stimulates the turnover rate of the QT-interval (S(D)). Solid arrows are production and loss processes, and dashed arrows are control actions.
The pharmacodynamic model is based on the approach presented in Guyton et al. (1972) and Snelder et al. (2014) consisting of three coupled turnover models describing HR, SV, and TPR, respectively. Each turnover model contains its individual turnover rate (kin) and fractional turnover rate (kout). The coupling of individual processes occurs via MAP. The time course of MAP is modeled as a nonlinear function according to eq. 2.
(2a)
(2b)CO denotes cardiac output, HRSV is a scaling constant, and BSLHR is the baseline value of heart rate. The term “baseline value” refers to the equilibrium value of a biomarker response (here HR) when no influences from handling, circadian variations, or drugs are present. This means that the baseline values are constant over time. MAP is coupled via feedback to an inhibitory action on HR, SV, and TPR. The pharmacodynamic model also includes circadian rhythm (CR) modeled as a 24-hour sinus function by means of eqs. 3 and 4:
(3)
(4)with ampHR, horHR, ampTPR, and horTPR representing amplitude and phase of HR and TPR. Handling effects (HDs) from the dosing procedures are modeled as an exponential process for HR
(5)and TPR
(6)with P being the amplitude, kHD being the time constant for the exponential decay, and t0 being handling time. CR and HD are assumed to affect the turnover rate (kin) of HR and TPR only, with no direct effect on SV. Drug intervention on HR, TPR, or QT is either stimulatory action (S(C)) given drug exposure (C)
(7)or inhibitory action (I(C)) given drug exposure (C) and metabolite concentration (Cm) (only present for pimobendan).
(8)For stimulatory effects on SV (from pimobendan), a linear model with a single slope parameter (SL) was used:
(9)Table 2 and eqs. 7 and 8 quantify each drug action by means of its own set of drug parameters: Emax/Imax, EC50/IC50, and γ.
The stimulatory (S) or inhibitory (I) action of each compound
The turnover rate of HR is expressed as
(10)in which turnover rate of HR is modeled by baseline turnover (kin,HR), circadian rhythm (CRHR), feedback (FB) from MAP, handling (HDHR), and fractional turnover rate (kout,HR). The drug effects (EFFHR) can be either stimulatory (EFFHR = S(C) eqs. 7) or inhibitory (EFFHR = I(C), eq. 8), with corresponding definitions for EFFTPR and EFFQT. EFFSV follows eq. 9. The turnover of the SV is expressed as
(11)The turnover rate of SV has a similar structure to that of HR but is assumed not to be affected directly by circadian rhythm or handling. The turnover of the TPR is expressed as
(12)which is assumed to have the same structure as that of HR. Initial conditions for eqs. 10–12 as well as eq. 18 below were determined numerically as described in the section Pharmacodynamic Model Parameters.
The FB from MAP is modeled as dependent on the individual baseline value (BSLMAP) for MAP, according to
(13)in which FB0 regulates the amount of feedback from MAP, and BSL0,MAP and FB0,MAP are empirically derived constants regulating the shape of the feedback relationship.
Turnover rates (kin) for the biomarker responses of HR, SV, and TPR were expressed as a function of their baseline values and the baseline value of MAP together with their fractional turnover rates (kout) and the feedback parameter (FB):
(14)The turnover rate of the SV was described according to
(15)The turnover rate of the TPR was described according to
(16)in which BSLHR, BSLSV, BSLTPR, and BSLMAP denote baseline values for HR, SV, TPR, and MAP, respectively. Here, we also require that
(17a)
(17b)The consistency of model equations with the definitions of the baseline values is further analyzed in Supplemental Section 6.
Modeling of QT Interval
The model was extended by adding the QT dynamics represented by a fourth turnover model. A feedback term was introduced to reflect the regulation of QT by means of HR, which is expressed according to
(18)Here, QTHR is a constant regulating the degree of QT adaptation to HR. The formulation of the handling effects, HDQT, follows those presented in eqs. 5 and 6:
(19)The turnover rate of the QT-interval was expressed as
(20)in which BSLQT and BSLHR denote baseline values for QT and HR, respectively.
Modeling of Concentration- and Biomarker Response–Time Data
Exploratory Analysis of High-Resolution Data.
An initial exploratory analysis of high-resolution and low-resolution biomarker-response data of HR, MAP, and QT-interval from the vehicle-control group included visual inspection. Biomarker responses displayed a slowly oscillating baseline response, which was repeatedly upwardly stimulated. A sinusoidal model (see Supplemental Section 1) was manually fitted to the slowly oscillating baseline response (for initial parameter estimates) until acceptable consistency was obtained between experimental and model-simulated data. Residuals between low-resolution data and sinus model predictions were computed, resulting in one set of residuals for HR and one set for MAP. These residuals were used to compute an initial distribution for the residual model in the NLME parameter estimation, as described below.
Drug Exposure Model Parameters.
Deviations between experimental and model predictions were assumed to follow a Gaussian distribution for all exposure data. A maximum likelihood objective function was used to produce a single estimate per compound of the pharmacokinetic parameters in eq. 1, a–c.
Pharmacodynamic Model Parameters.
The parameters of the pharmacodynamic model were estimated using an NLME approach, allowing simultaneous regression of all individual biomarker-response time courses (i.e., HR, MAP, and QT). Note that the stroke volume was not among the biomarkers measured in this work. In the cases in which the model was regressed to effects on stroke volume (pimobendan), it was done based on the secondary effects on the measured biomarkers. This technique allows the estimation of interindividual variability directly from data. Specifically, the analysis was done using a custom implementation of the expectation-maximization algorithm (Kuhn and Lavielle, 2005) in MATLAB (version R2019b, The Mathworks, Inc., Natick, Mass.) to maximize the likelihood, L, of the system and drug-specific parameters (
) given the measured data, 
:
(21)Here, N denotes the total number of individuals, and ηi denotes the vector of individual random-effects model parameters for individual (i). The distribution of the individual random effects was assumed to be Gaussian. For increased stability of the estimation algorithm, 1000 sampled parameter vectors for each individual were used in each iteration. This also obviated the need for smoothing the estimated parameter values at the cost of a higher computational demand. Convergence was considered to be attained when none of the parameters diverged more than 1% from their average across 10 iterations of the algorithm.
The baseline parameters BSLHR, BSLMAP, and BSLQT were assumed to vary across individuals. No interindividual variability was added to drug parameters Emax, EC50, and γ of the modelled biomarker responses (HR, SV, TPR, and QT). A Gaussian prior with a mean of 0.126 and a S.D. of 0.013 was used to constrain kout,SV to values around previously reported values for rat (Snelder et al., 2014) since this was found to improve estimation stability. After a change in any parameter value, new initial conditions were computed by initializing each modeled biomarker to its BSL and simulating the model over a long time span (336 hours), allowing it to reach equilibrium.
Equilibrium Concentration-Response Relationships.
The estimated pharmacokinetic/pharmacodynamic model was used to predict the equilibrium concentration–biomarker response relationship of each compound. Here, equilibrium refers to a state in which the unbound plasma concentration is constant and wherein no circadian rhythm is present. Since the model is nonlinear, the biomarker response cannot be expressed explicitly in terms of drug exposure. Instead, steady-state data were generated by means of model simulations (t = 96 hours) with a logarithmically spaced exposure range. Parameters for handling (PHR, PTPR, PQT), diurnal variations (ampHR, ampTPR), and residual errors were set to zero. Standard target-effect models for sigmoidal Emax and Imax models with baseline [Table 3.1, p.221 in Gabrielsson and Weiner (2016)], according to eqs. 22a or 22b below, were then regressed to the simulated concentration–biomarker response data. In other words,
(22a)
(22b)depending on whether the observed effect was of stimulatory or inhibitory action, respectively.
From the estimated drug parameters EC50, Emax, and γ at the equilibrium state, an effective target concentration (CT) was calculated using the standard target concentration model for a sigmoidal Emax model with baseline wherein CT = EC50 ((ET – E0)/(E0 + Emax – ET))1/γ, which for a 20% change in response from baseline (i.e., ET = E0 + 0.2 Emax) results in
(23a)and likewise, in the case of inhibitory action with ET = E0 – 0.2 Imax
(23b)This analysis was not applicable to the effects on SV from pimobendan since these were modeled as linear according to eq. 9.
Results
Exploratory Analysis of High-Resolution Data.
High-resolution biomarker-response data (obtained from ECG and blood pressure measured at 1000 Hz, downsampled to 500 Hz) from four rats were analyzed in step 1 to understand baseline behavior and residual analysis. An example of high-resolution response of HR and MAP for a single animal obtained during a 24-hour vehicle-control study is shown in Fig. 4.
Upper panel: high-resolution response-time data of heart rate during a 24-hour baseline cycle. Lower panel: high-resolution response-time data of blood pressure (MAP) during a 24-hour baseline cycle. Dashed vertical black lines indicate start and stop of the 12-hour light/dark cycle. Red arrow shows time of oral (saline) dosing.
HR and MAP responses randomly deviated toward higher values and were superimposed on a diurnal baseline wave (example results for one of the animals are shown in Fig. 5). This is particularly evident for HR (upper panel, Fig. 4) and MAP (lower panel, Fig. 4). Light and dark periods differed clearly in duration, frequency, and amplitude of deviations. A sinusoidal model was fitted to the diurnal baseline wave obtained from high-resolution data sets. The residuals between the model-predicted and experimental data were obtained for low-resolution data. The high- and low-resolution responses and model regressions are compared in the upper panel of Fig. 5.
Upper row: observed heart rate and blood pressure (MAP) vs. clock time of high- and low-resolution data (gray and orange lines, respectively) with a sinusoidal model fit to the resting state (black solid line). The light and dark periods are separated by vertical dashed lines. Bottom row: observed blood pressure (MAP) vs. heart rate residuals (solid symbols) between model and the low-resolution data. The solid (blue) lines represent fitted log-normal distributions of light (left) and dark (right) periods.
To capture both the asymmetry of HR and MAP around their respective baselines and the correlation between the two signals the bivariate log-normal distribution was chosen to describe the residual errors. Because of the systematic differences between light- and dark-hour response-time data, two separate residual error models were used to describe light and dark conditions (depicted in Fig. 5, bottom panel). Parameters for the residual error model are given in Supplemental Section 2. For the QT-interval response, a normal distribution was used since no correlation to either HR or MAP (unpublished data) was found (see Supplemental Fig. 1). These residual error models were the main outcome of the exploratory analysis and were used in the pharmacodynamics analysis.
Exposure Time Courses of Test Compounds.
The two-compartment disposition model (Fig. 2; eqs. 1, a–c) was simultaneously fitted to intravenous and oral data shown in Fig. 6 for each compound. Final parameter estimates of disposition (CL, CLd, and volumes of the central and peripheral compartments) and absorption parameters (Ka and bioavailability (F)) are given in Table 3.
Observed total concentration-time data after intravenous (solid triangles) and oral (solid squares) dosing of ivabradine, sildenafil, dofetilide, and pimobendan. Solid lines are model-fitted concentration time courses driving the pharmacodynamic model. Doses of ivabradine, sildenafil, dofetilide, and pimobendan were 2.13, 21.1, 2.26, and 8.97 µmolkg−1 (oral) and 1, 1, 10, and 1 µmolkg−1 (intravenous), respectively. The horizontal dashed red lines represent unbound potencies in Table 4 and are shown as a comparison between exposure data and the predicted drug parameter. The red, simulated plasma concentration-time curves are model simulations of the highest oral dose of each compound in the pharmacodynamic assessment of biomarker data.
Final pharmacokinetic parameter estimates and their RSE%
The goodness of fit was assessed by means of residual analysis. Note that several of the parameters are poorly estimated (relative standard error (RSE) >100%), which is not commonly accepted for correct assessment of the pharmacokinetic properties. However, the goal was to correctly capture both high- and low-exposure data to have a pharmacokinetic model that predicts exposure data well and can serve as a “driver” of biomarker-response data. Concentration time courses of experimental and model-predicted data are shown in Fig. 6. The oral model was then used to predict the time courses of unbound plasma concentrations for the actual dosing regimens used in step 2 of pharmacodynamic assessment of data. A common concentration time course was used for all animals in each drug study.
Pharmacodynamic Model of Biomarker Responses.
The agreement between observed and model-predicted data was assessed using predictive plots presented in Figs. 7–10. The data shown in the predictive plots were computed from model simulations based on 1000 parameter sets sampled from the estimated population distribution. Then, 1000 predictions of the residuals each of HR, MAP, and QT were sampled and added to the model output. Finally, for each time point, the 1000 simulated responses were ordered, and the 25 highest and lowest values were removed to construct a 95% prediction interval. Regressed exposure data showed high consistency between measured data and model output, with the majority of data points well within their prediction intervals. However, the model appears to systematically overpredict the drug blood pressure effect of the highest pimobendan dose.
Experimental data (gray dots), model-predicted avg. (solid line), and 95% prediction interval (dashed lines) response of heart rate and blood pressure (MAP) after a single oral dose of ivabradine at 0, 3, 10, and 30 mgkg−1 corresponding to 0, 6.40, 21.3, and 64.0 µmolkg−1. The predicted vehicle-control response (gray line) is superimposed on experimental data.
Experimental data (gray dots), model-predicted avg. (solid line), and 95% prediction interval (dashed lines) response of heart rate and blood pressure (MAP) after a single oral dose of sildenafil at 0, 3, 10, and 30 mgkg−1 corresponding to 0, 6.32, 21.1, and 63.2 µmolkg−1. The predicted vehicle-control response (gray line) is superimposed on experimental data.
Experimental data (gray dots), model-predicted avg. (solid line), and 95% prediction interval (dashed lines) response of QT-interval after a single oral dose of dofetilide at 0, 3, 10, and 30 mgkg−1 corresponding to 0, 6.79, 22.6, and 67.9 µmolkg−1. The predicted vehicle-control response (gray line) is superimposed on experimental data.
Experimental data (gray dots), model-predicted avg. (solid line), and 95% prediction interval (dashed lines) response of heart rate and blood pressure (MAP) after a single oral dose of pimobendan at 0, 3, 10, and 30 mgkg−1 corresponding to 0, 8.96, 29.9, and 89.7 µmolkg−1. The predicted vehicle-control response (gray line) is superimposed on experimental data.
Parameter estimates were obtained from low-resolution data, utilizing residual distributions for light and dark periods from the exploratory analysis of the high-resolution data as initial estimates. Estimated pharmacodynamic parameters and their RSE% are presented in Table 4. Additional fixed-model parameters are listed in Supplemental Table 1.
Estimated systems and drug parameters in rats and their precision (RSE%)
All IC50/EC50 refer to unbound conc.
For the system’s parameters, estimated baseline values for heart rate (BSLHR), mean arterial pressure (BSLMAP), QT-interval (BSLQT), the parameter regulating baroreflex feedback (FB0), and handling (PHR and PTPR) were estimated with good precision. The available data exhibited a very weak relationship between HR and QT (see Supplemental Fig. 1), which led us to exclude the HRQT from the estimation by fixing it to zero. The interindividual variability was low (5%). Drug-specific parameters, such as Emax, EC50, and γ, had generally high precision except for the pimobendan action on stroke volume (RSE% >100%) and dofetilide IC50 (RSE% = 93%). The potencies of ivabradine (HR) and pimobendan (SV) were estimated to be 475 and 604 nM, respectively, which is beyond the observed exposure range but still within published data (Kitzen et al., 1988; Du et al., 2004). The potencies of dofetilide (QT) and sildenafil (TPR) were estimated to be 50.6 and 4.01 nM, respectively, which is within the observed exposure range and consistent with published data (Mounsey and DiMarco, 2000; Gresser and Gleiter, 2002). The γ parameter was relatively low for all compounds and biomarker responses and particularly for dofetilide and QT effect, suggesting a shallow equilibrium concentration–biomarker response relationship around its potency value (Table 5).
Estimated equilibrium parameters and predicted EC20 values in rats and free conc. Cu in human
t1/2,kout from Table 4 are included for convenience. Cu denotes clinical plasma conc. in humans. All conc. are unbound.
Equilibrium Concentration–Biomarker Response Analyses.
The equilibrium concentration–biomarker response (heart rate, blood pressure, and QT-interval) relationships were assessed for each test compound (Figs. 11–13). Supplemental Section 5 shows how data were generated.
Model-predicted equilibrium concentration–heart rate response (solid red line), population variability (black dotted lines), and 20% response in heart rate (black dashed line) for ivabradine, sildenafil, dofetilide, and pimobendan, respectively.
Model-predicted equilibrium concentration–blood pressure (MAP) response (solid red line), population variability (black dotted lines), and 20% change in blood pressure (black dashed line) for ivabradine, sildenafil, dofetilide, and pimobendan, respectively.
Model-predicted equilibrium concentration–QT response (solid red line), population variability (black dotted lines), and 20% change in QT response (black dashed line) for ivabradine, sildenafil, dofetilide, and pimobendan, respectively.
Emax/Imax, EC50/IC50, and γ of the equilibrium concentration–response relationships were estimated from model-generated data in Table 4 (Figs. 11–13; Table 5) by regressing eq. 22a or eq. 22b to the simulated data. In the regression, the values of E0 in these equations correspond to the baseline values BSLHR, BSLQT, and BSLMAP listed in Table 4. The clinical free concentration (Cu) ranges are included for comparisons.
Discussion
This analysis covers exposure- and response-time data from HR and MAP data after acute oral dosing of ivabradine, sildenafil, and pimobendan in Han-Wistar rats. Data on QT-interval duration were also obtained from dofetilide. The analysis allowed system parameters (turnover rates, etc.) to be shared across test compounds and studies, whereas drug parameters (such as EC50) were compound-specific and shared across dosing regimens. Finally, the equilibrium concentration–response relationships were established and visualized.
Exploratory Analysis of High-Resolution Biomarker-Response Data.
The exploratory data analysis covered both high- and low-resolution biomarker-response data, revealing two structural features. One was a slowly oscillating baseline, which was assumed to originate from intrinsic circadian variations in biomarker response. The second feature was that fluctuations from baseline increased both in amplitude and duration during dark periods, which are believed to originate from acute baroreflex resetting due to increased animal activity (Potts and Mitchell, 1998; Dampney, 2017). Data suggested that a model of residual errors would be asymmetric and correlated, which led to selection of a translated log-normal distribution model with dark/light-specific parameters.
Exposure Time Courses of Test Compounds.
The exposure time courses of ivabradine (Zhang et al., 2016), sildenafil (Sawatdee et al., 2018), dofetilide (Smith et al., 1992), and pimobendan (Asakura et al., 1993) were previously studied in Sprague-Dawley rats (200–300 g) and showed slightly different pharmacokinetic properties compared with the larger Han-Wistar rats (450–600 g) used here (Fig. 6). This highlights the necessity of exposure time courses obtained in the actual animal strain for model development. The primary goal was to apply an exposure model that accurately captures both high– and low–exposure time data for each test compound.
Pharmacodynamic Model.
A series of pharmacodynamic models of cardiovascular effects in rodents have been developed for BP- (Hao et al., 2007; Bertera et al., 2012; Kiriyama et al., 2016), HR- (van Steeg et al., 2007; Bertera et al., 2012; Kiriyama et al., 2016), and ECG-based biomarkers (Bol et al., 1997; Ohtani et al., 2000; Kiriyama et al., 2016). The majority of these neglect the intertwined dependencies across cardiovascular measures, making comparison of results difficult. For instance, γ parameters in Tables 4 and 5 were generally lower than previously reported values in human (Chu et al., 1999; Duffull et al., 2000; Mehrotra et al., 2007; Gotta et al., 2016), which were generally reported to be ≥1, but it is unclear whether this stems from differences in species or model. A few studies (Francheteau et al., 1993; Upton and Ludbrook, 2005; Snelder et al., 2014; Kamendi et al., 2016) have applied an integrative approach, including inter-relationships between TPR, SV, and MAP.
Estimated baseline values of MAP, HR, and QT and their variabilities are consistent with previously published data (Howgate, 2013; Snelder et al., 2014). The estimated feedback coefficient FB0 shows a more effective regulation in Han-Wistar rats compared with Sprague-Dawley rats (0.0053 vs. 0.0029). Predicted half-lives based on estimated fractional turnover rates (kout) were consistent with published data: 3.6 minutes of HR, 5.5 hours of SV, and 12 minutes of TPR. This suggests more rapid dynamics in HR, QT, and SV and slower dynamics in TPR. All system parameters were estimated with good precision (1.1%–21% RSE). Although HRQT was excluded from estimation because of the weak HR-QT relationship in rodents, it was included in the model structure for completeness and in anticipation of studies in nonrodent species.
Table 4 shows estimated parameters related to their specific action on HR, SV, TPR, or QT (eqs. 10–12, and 18). In contrast, Table 5 shows estimates obtained from the predicted equilibrium concentration–response relationship based on simulated steady-state data. Thus, values in Table 5 take into account the systemic interaction between the biomarkers via MAP (eqs. 2 and 10–12), whereas values in Table 4 do not.
Ivabradine is a selective blocker of the funny current that specifically reduces heart rate (Tardif et al., 2009). However, with high doses, the drastic negative chronotropic effects also affect the cardiac output, lowering the arterial blood pressure (Fig. 7). Because of the indirect effect on blood pressure, only drug-specific parameters for heart rate are presented in Table 4. However, the steady-state analysis presented in Table 5 provides potencies for both heart rate and blood pressure. The unbound plasma concentration in Table 5 shows an IC20 (HR) close to clinical peak (unbound) concentration, whereas the corresponding value for IC20 (MAP) demonstrates a 300-fold margin. Although estimated steady-state IC50 values for heart rate are higher than values reported for human (Duffull et al., 2000), they are 100-fold lower than corresponding values for MAP, being consistent with the observation that ivabradine does not give a significant decrease in blood pressure at therapeutic plasma concentrations (Deedwania, 2013; Koruth et al., 2017). Despite a pronounced reduction in heart rate at high doses, no relevant change in QT-interval duration was observed. A weak relationship between QT and RR intervals has been previously reported in mice (Roussel et al., 2016) and more recently in unanesthetized rats (Mulla et al., 2018), suggesting that the influence of HR on QT is marginal in rodents.
Sildenafil promotes peripheral vasodilation. Simulations (Fig. 8) showed a marked increase in HR and a minor reduction in MAP, consistent with reported results at therapeutic concentrations (Vardi et al., 2002). For sildenafil, EC20 (HR) and IC50 (MAP) were reached by a 5-fold increase from the lowest Cu, which was close to the upper limit of the therapeutic interval (Table 5).
Dofetilide is a highly potent blocker of the rapid component of the IKr current carried by hERG channel in humans, thereby increasing the action-potential duration and QT-interval (Trudeau et al., 1995). In this study, dofetilide induced a pronounced prolongation of the QT-interval duration, with steady-state values for both EC20 and EC50 well below the free peak concentration Cmax of the lowest dose. The rat is controversial as model for assessing QT prolongation because of its low expression or weak functionality of hERG-like channels in the ventricles (Rees and Curtis, 1996; McDermott et al., 2002). However, expression of ERG in rat heart tissue has been indicated using RNase protection assays and confirmed by blocking the tail current by E-4031 and dofetilide in rat myocytes (Wymore et al., 1997). Likewise, compounds with different modes of action (e.g., erythromycin) have been observed to affect QT prolongation in rodents (Ohtani et al., 2000; McDermott et al., 2002). These findings suggest that the contribution of the IKr component to the ventricular repolarization in rodents is present but weak, making the detection of QT prolongation in vivo challenging. The Hill coefficient of the dofetilide exposure–QT response relationship was low, suggesting a shallow equilibrium concentration–response relationship but, more importantly, an extended duration of QT response whenever the effect has been established in spite of rapidly declining plasma exposure. On the other hand, below critical exposure at the target site, rapid plasma concentration fluctuations may not be deleterious because target exposure will always lag after the plasma exposure by a half-life of about 4 minutes (Table 5). The time to establish equilibrium between plasma and target-site exposure will be 3 to 4 half-lives (10–15 minutes).
Pimobendan is a calcium sensitizer and a selective PDE3 inhibitor with positive inotropic and vasodilator effects primarily used in veterinary medicine (Kitzen et al., 1988). Estimated values for the slope parameter SLSV (inotropic effects) and IC50,TPR (vasodilator effects) show a higher potency for vasodilation compared with stroke volume, which are negligible within the investigated exposure range (Table 4). This suggests that the model impacts the vasodilator effect as a cause of the observed changes in blood pressure and heart rate. However, Fig. 10 and Supplemental Fig. 9 show that the model has difficulties adapting to high-dose effects on blood pressure. The EC20/IC20 estimates correspond to a 9-fold increase from the lowest Cu. Pimobendan has been reported to lower MAP and increase HR in several animal species (Kitzen et al., 1988) but to a lesser extent in humans (Kubo et al., 1992; Chu et al., 1999).
Although the model is able to mimic both vehicle-control and drug responses, some systematic deviations between measured and modeled percentiles can be observed in Supplemental Figs. 3–9 (Supplemental Material). This might indicate a need for further model refinement and that extrapolating model predictions should be done with caution.
New Features of the Proposed Model.
The current approach in safety pharmacology studies is to test compounds as a single dose at two to three dose levels defined to ensure that plasma exposure is at least 3- to 30-fold greater than the estimated therapeutic peak exposure. We extended previously published models to also incorporate turnover of QT response (Fig. 2). This allows quantification of pivotal safety measures, such as the margin between therapeutic and extended exposure with a risk of adverse effects, exemplified by the equilibrium EC20 and EC50. The rationale for including the QT-interval achieved in rats is primarily of translational safety reasons. This study constitutes a first evaluation of the model on rat data, but further applications of the proposed model on data from other species are ongoing.
The exploratory analysis revealed a more rational description of shifting biomarker responses during light and dark conditions, allowing the model to naturally replicate previously observed asymmetric variations in HR and MAP. The integrated model was used to predict the concentration-response relationships of the four reference compounds and to quantify readouts pivotal for safety assessment, such as onset, intensity, and duration of response as well as the equilibrium exposure–response relationships. Here, EC20 was used as an example safety measure, but any appropriate measure can be used within the presented methodology.
The proposed model may be applied to highlight the plasma exposure associated with potential adverse effects within and beyond therapeutic exposure and predictions of safety margins without necessarily additional preclinical studies. This analysis demonstrates the utility of a model-based approach that integrates data from different studies and compounds for refined preclinical assessment of safety margins.
Acknowledgments
The authors are grateful to Werner Mayer, who performed the in vivo experiments in rats, and to Hermann Rapp, who provided valuable assistance with PK-related issues.
Authorship Contributions
Participated in research design: Scheuerer, Martel, Pairet.
Conducted experiments: Scheuerer, Pairet.
Performed data analysis: Wallman, Jirstrand, Gabrielsson.
Wrote or contributed to the writing of the manuscript: Wallman, Scheuerer, Martel, Pairet, Jirstrand, Gabrielsson.
Footnotes
- Received September 21, 2020.
- Accepted February 22, 2021.
This work received no external funding.
Conflict of Interest Statement: No author has an actual or perceived conflict of interest with the contents of this article.
↵
This article has supplemental material available at jpet.aspetjournals.org.
Abbreviations
- BP
- blood pressure
- BSL
- baseline
- CL
- clearance
- CLd
- intercompartmental distribution
- Cu
- clinical free concentration
- DSI
- Data Science International
- EC20
- 20% effective concentration
- ERG
- ether-a-go-go–related gene
- F
- bioavailability
- FB
- feedback
- hERG
- human ERG
- HR
- heart rate
- IC20
- 20% inhibitory concentration
- IKr
- rapid delayed rectifier current
- Ka
- absorption rate constant
- kin
- turnover rate
- kout
- fractional turnover rate
- MAP
- mean arterial BP
- NLME
- nonlinear mixed-effects
- PDE
- phosphodiesterase
- PK
- pharmacokinetic
- RSE
- relative standard error
- SV
- stroke volume
- TPR
- total peripheral resistance
- Copyright © 2021 by The Author(s)
This is an open access article distributed under the CC BY-NC Attribution 4.0 International license.
References
In this issue
Improved Assessment of Cardiovascular Safety Data
