Joint feedback analysis modeling of nonesterified fatty acids in obese Zucker rats and normal Sprague-Dawley rats after different routes of administration of nicotinic acid.

Data were pooled from several studies on nicotinic acid (NiAc) intervention of fatty acid turnover in normal Sprague-Dawley and obese Zucker rats in order to perform a joint PKPD of data from more than 100 normal Sprague-Dawley and obese Zucker rats, exposed to several administration routes and rates. To describe the difference in pharmacodynamic parameters between obese and normal rats, we modified a previously published nonlinear mixed effects model describing tolerance and oscillatory rebound effects of NiAc on nonesterified fatty acids plasma concentrations. An important conclusion is that planning of experiments and dose scheduling cannot rely on pilot studies on normal animals alone. The obese rats have a less-pronounced concentration-response relationship and need higher doses to exhibit desired response. The relative level of fatty acid rebound after cessation of NiAc administration was also quantified in the two rat populations. Building joint normal-disease models with scaling parameter(s) to characterize the "degree of disease" can be a useful tool when designing informative experiments on diseased animals, particularly in the preclinical screen. Data were analyzed using nonlinear mixed effects modeling, for the optimization, we used an improved method for calculating the gradient than the usually adopted finite difference approximation.

Introduction

Nicotinic acid (NiAc) is a lipid lowering agent that inhibits lipolysis in adipose tissue by activating the GPR109A receptor, resulting in a pronounced decrease in plasma nonesterified fatty acids (NEFAs) concentrations.1 Current understanding of the mechanism is that activation of the G protein-coupled receptor GPR109A by NiAc inhibits adenylyl cyclase activity, leading to decreased formation of cyclic adenosine monophosphate (cAMP) from adenosine triphosphate. cAMP regulates lipolysis in adipocytes by activating protein kinase A that phosphorylates hormone-sensitive lipase (HSL). NiAc thus reduces the hydrolysis of triglycerides (TG) into NEFAs and glycerol, which is catalyzed by HSL.2 A feedback model describing NiAc-induced changes in NEFA plasma concentrations has previously been published.3–5 The characteristics described by the model are tolerance and oscillatory rebound effects, and depending on the parameter values, the model can be used for either Sprague–Dawley rats or obese Zucker rats.

Obese Zucker rats display insulin resistance and obesity, and are frequently used as a disease model. Obesity is frequently associated with insulin resistance and known to influence the distribution and clearance of compounds.6–10 This suggests that insulin resistance, in an animal or a patient, can alter the pharmacokinetics (PK) and pharmacodynamics (PD) of the NiAc/NEFA system with changes in onset, duration and intensity of the drug effects. As a normal model, Sprague–Dawley rats are used. Their lipoprotein and lipid patterns are similar to those of the lean Zucker rat and therefore can serve as a normal model for comparison purposes.11,12

Analyzing all data collected during a preclinical screen allows us to better quantify the difference in PD between normal Sprague–Dawley rats and obese Zucker rats. A full analysis of all data available also allows us to better understand how drug response in normal animals translates to the diseased animals. This is particularly valuable for decisions made in the preclinical phase, which are often based on selected data from a normal animal model. We therefore sought to determine how disease affects the PD and to quantify this effect. This report describes such a joint analysis that utilized data from 95 Sprague–Dawley rats and 19 obese Zucker rats simultaneously, to estimate the relative change in pivotal parameters across normal and diseased animals. Data were analyzed using nonlinear mixed effects modeling in which, for the optimization, we used a more exact method for calculating the gradient than the usually adopted finite difference approximation.13

We also compare our joint analysis to population PK/PD analysis carried out separately for normal and obese rats to investigate which PD parameters are affected by the disease. The combined analysis uses all available data simultaneously to maximize the information content from which inference about a pathophysiological system can be made. In a more general setting, this kind of modeling can be used to combine knowledge gained from previous, less information-rich studies. Accumulating information into models over the life-span of a discovery project, can in this way reduce the number of future measurements or studies but still keep the quantitative quality of results.

Materials and Methods

The following material is mainly adapted from previous publications3–5 but is repeated here for the convenience of the reader.

Chemicals

Nicotinic acid (pyridine-3-carboxylic acid) was obtained from Sigma–Aldrich (St. Louis, Missouri) and was dissolved in 0.9% NaCl. All solvents were of analytical grade and the water used in the experiments was obtained from a water purification system (Elgastat Maxima, ELGA, Lane End, UK).

Animals and Surgical Procedures

Male obese Zucker (fa/fa) and normal male Sprague–Dawley rats were purchased from Harlan Laboratories Nederlands B.V. (Horst, the Netherlands) at 7 and 11 weeks of age, respectively, and used at 16 weeks of age. The animals were housed in groups of 5–6 with free access to standard rodent chow (R3; Laktamin AB, Stockholm, Sweden) and tap water. They were kept in climate-controlled facilities at a room temperature of 20°C–22°C and relative humidity of 40%–60% under a 12:12-h light–dark cycle. The study was approved by the Ethics Committee for Animal Experiments, Gothenburg, Sweden (EA 100868). Surgery was performed under isoflurane (ForeneQR; Abbott Scandinavia AB, Solna, Sweden) anesthesia and body temperature was maintained at 37°C using a thermoregulated heating pad. Catheters were implanted in the left carotid artery for blood sampling and in the right external jugular vein for drug administration, as previously described.14 Before cannulation, catheters (IntramedicQR, PE50; Becton, Dickinson and Company, Franklin Lakes, New Jersey) were filled with sterile sodium-citrate solution (20.6 mM sodium-citrate in sterile saline; Pharmaceutical and Analytical R&D, AstraZeneca, Mölndal, Sweden) to prevent clotting. After cannulation, the catheters were exteriorized at the nape of the neck and sealed. After surgery, the rats were housed individually and allowed 5 days to recover before the experiments began.

Experimental Design

The animals were fasted for 14 h before dosing and throughout the experiment to minimize the fluctuations in NEFA caused by food intake. They had free access to drinking water during the length of fast. On the day of experimentation, they were weighed, moved to clean cages, and the venous catheters were connected to infusion pumps (CMA 100; Carnegie Medicin AB, Stockholm, Sweden). Following a 30 min adaptation period, two consecutive arterial blood samples were collected 15 and 5 min before drug administration to determine predose baseline NEFA and NiAc concentrations. Normal Sprague–Dawley rats (weighing 220–367 g) were assigned to 12 groups. Groups 1–8 received an intravenous constant rate infusion for either 30 or 300 min. Four of the eight groups received vehicle (0.9% NaCl, n = 10), or 1 (n = 4), 5 (n = 8) or 20 (n = 9) μmol kg−1 NiAc over 30 min. The remaining four groups received vehicle (n = 8), or 5 (n = 9), 10 (n = 8) or 51 (n = 7) μmol kg−1 NiAc over 300 min. Groups 9–12 received oral doses by gavage with vehicle or 24.4, 81.2, or 812 μmol kg−1 NiAc (n = 6 per group). The concentrations of the dosing solutions were adjusted to give infusion volume flow rates in the range of 0.4–22 μL min−1 and oral dosing volume in the range of 1.4–1.6 mL, based on body weight. The dosing solutions were prepared within 30 min of administration by dissolving an appropriate amount of NiAc in saline solution. Among groups 9–12, two groups (n = 5 in each) received a constant infusion of NiAc of 5 μmol kg−1 for 30 min, followed by a stepwise decrease in infusion rate every 10 min. The last of the two groups was given a 5 μmol kg−1 for 30 min at 210 min. Control groups had the same schemes (n = 1 in each) but received vehicle. Obese Zucker rats (weighing 473–547 g) were assigned to four groups of which received an intravenous constant rate infusion for either 30 or 300 min. Two groups received vehicle (0.9% NaCl, n = 2) or 20 (n = 8) μmol kg−1 NiAc over 30 min. The other two groups received vehicle (n = 2) or 51 (n = 7) μmol kg−1 NiAc over 300 min.

Multiple arterial blood samples were drawn, 11–15 per rat, for both the 30 and 300 min infusion experiments for analysis of NiAc and NEFA plasma concentrations. The total blood volume removed did not exceed 1.5 mL, and was replaced with an equal volume of sterile sodium-citrate solution to maintain a constant circulatory volume. The control rats received the same volume of infusion solution (vehicle) as the NiAc groups, and all animals were subjected to similar sampling procedures. The blood samples (120 μL each) were collected in EDTA coated polyethylene tubes and kept on ice until centrifuged (10,000g, 5 min, 4°C). The plasma was stored at −20°C pending analysis. The start of infusion was taken as time zero (0 min).

Analytical Assays

Analysis and quantification of NiAc in plasma were performed using liquid chromatography—mass spectrometry. The HPLC system was an Agilent 1100 Series (Hewlett-Packard GmbH, Walbronn, Germany) coupled to an HTC PAL autosampler (CTC Analytics AG, Zwingen, Germany). Plasma samples (50 μL per sample) were precipitated with cold acetonitrile containing 0.2% formic acid (150 μL per sample). After vortex mixing and centrifugation at 4°C (4000g, 20 min), an aliquot of 100 μL of the supernatant was used for analysis. The mobile phase consisted of (A) 2% acetonitrile and 0.2% formic acid in water, and (B) 0.2% formic acid in acetonitrile. Separation was performed on a 50 × 2.1 mm Biobasic AX column with 5 lm particles (Thermo Hypersil-Keystone, Runcorn, Cheshire, UK) with a gradient of 95%–20% B over 1 min, held at 20% B for 1.5 min, and returned to initial conditions in one step. The HPLC system was connected to a Sciex API 4000 quadrupole mass spectrometer with a positive electrospray ionization interface (Applied Biosystems, Ontario, Canada) and the mass transition was 124.0 > 80.2. Data acquisition and data evaluation were performed using Analyst 1.4.1 (Applied Biosystems). The method showed linearity over a concentration range of 0.001–28 μmol L−1. The lower limit of quantification was 1 nmol L−1 applying a sample volume of 50 μL plasma. Plasma NEFA was analyzed using an enzymatic colorimetric method (Wako Chemicals GmbH, Neuss, Germany) adapted to a 96-well format.5

Pharmacokinetic Model of NiAc Disposition

Pharmacokinetic models of different complexity have previously been evaluated3–5 for both normal and obese animals. As this study primarily focus on quantifying differences in PD across different animal populations, the PK model used to drive the PD model was adapted without alterations from a previous study.3,15 A schematic picture of the chosen models is shown in Figure 1.

Figure 1

(a) Schematic model of the absorption and NiAc disposition. Ag is the amount of NiAc in the gut; Cp and Ct are the plasma and peripheral concentrations, respectively. Dosep denotes the oral dosing, Synt the endogenous synthesis of NiAc, and Inf the intravenous infusion. NiAc disposition in normal rats is described by a two compartment model with two nonlinear elimination terms. Translocation from the gastrointestinal tract is described by a linear and a nonlinear process. (b) Schematic model of NiAc disposition for obese Zucker rats. NiAc disposition is described in obese rats by a one compartment model with a nonlinear elimination term.

The PK model handles the different administration routes, such as constant infusion, stepwise infusion schemes, and oral dosing. For the oral dosing, disappearance of NiAc from the gastrointestinal tract was modeled by one linear and one capacity-limited elimination process operating in parallel according to

where Ag is the amount of drug in the gut. The rates A1 and A2 correspond to the linear and the nonlinear translocation of the drug from the gut, respectively. The parameter ka is a first order absorption constant, Vmax,g is the maximum absorption rate, and Km,g corresponds to the amount of drug in the gut for which the rate of the capacity-limited elimination process is 50% of Vmax,g. The disposition model of NiAc for normal animals is described by

where Cp and Ct denote the NiAc concentration in the central and peripheral compartments with volumes Vc and Vt, Inf the drug infusion rate of NiAc, Synt the endogenous synthesis rate, Vmax1 and Km1 the maximal rate and Michaelis–Menten parameter of the high affinity process, Vmax2 and Km2 the maximal rate and Michaelis–Menten parameter of the low affinity process, and Cld the intercompartmental clearance. The endogenous concentration of NiAc was estimated according to the solution of the Eqs. 2 and 3 at steady state. The disposition of NiAc in obese Zucker rats was modeled by a one compartment model with endogenous synthesis, Synt, of NiAc and capacity-limited elimination (see Fig. 1)

where Cp denotes the NiAc concentration in the central compartment, Vc the central volume of distribution, Inf the drug infusion rate, Synt the endogenous synthesis rate. The parameters Vmax1 and Km1 denote the maximal rate and the Michaelis–Menten parameter, respectively. In Table1, the previously estimated model parameter values are shown for both normal and obese rats.

Table 1

Fixed Effects Population Pharmacokinetic Parameter Estimates and Interindividual Variability (IIV) with Corresponding Relative Standard Errors (RSE%)

		Obese Zucker Rats		Normal Sprague–Dawley Rats
Parameter	Definitions	Estimate (RSE%)	IIV (RSE%)	Estimate (RSE%)	IIV (RSE%)
Vmax1 (μmol min−1 kg−1)	Max. velocity, pathway 1	1.59 (13.9)	21.4 (234)	0.0871 (22.8)	92.7 (27.5)
Km1(μmol L−1)	MMa constant, pathway 1	18.9 (21.5)	–	0.235 (29.2)	–
Vmax2 (μmol min−1 kg−1)	Max. velocity, pathway 2	–	–	7.09 (39.6)	29.1 (43.6)
Km2 (μmol L−1)	MMa constant, pathway 2	–	–	74.5 (43.4)	–
Vc (L kg−1)	Central volume	0.323 (12.4)	–	0.393 (5.29)	–
Vt (L kg−1)	Peripheral volume	–	–	0.172 (35.2)	–
Cld (L min−1 kg−1)	Intercompartmental distr.	–	–	0.000852 (27.8)	–
Synt (μmol min−1 kg−1)	Endogenous synthesis rate	0.00280 (10.1)	95.3 (115)	0.00355 (23.3)	109 (34.7)
ka (min−1)	First order absorption rate	–	–	0.00477 (33.5)	10.8 (93.1)
Vmax,g (μmol min−1kg−1)	Maximum absorption rate	–	–	2.96 (21.3)	10.8 (93.1)
Km,g (μmol kg−1)	Amt. in gut at half Vmax,g	–	–	20.5 (29.6)	–
σ1	Residual prop. Error	40.0 (26.3)	–	42.8 (5.16)	–
σ2	Residual add. Error	–	–	–	–

Previously published in Refs. 4–5.

Michaelis–Menten constant.

Feedback Model of NEFA

To describe the indirect response relationship between NiAc and NEFA plasma concentrations, the feedback model in Figure 2 is used. NiAc affects NEFA plasma concentrations by inhibition of hydrolysis of TG to NEFA and glycerol in adipocytes, thereby reducing the release of NEFA into plasma.16–18 A fraction of NEFA will remain in the circulation, unaffected by NiAc, representing the lower physiological limit of NEFA in plasma. This process is incorporated as a zero-order production term, kcap, in the model. Feedback in the system was modeled as a moderator distributed over a series of eight transit compartments where the moderator in the first compartment, M1, inhibited the formation of NEFA, which we denote by R (i.e., the build-up of NEFA concentrations). NEFA formation is modeled as

Figure 2

Schematic feedback model of NEFA production in normal and obese rats. The NEFA turnover is described as a feedback model with eight moderator transit compartments. The solid and dashed lines represent fluxes and control processes, respectively. See Table2 for an explanation of each model parameter.

where M1 is the moderator in compartment one, M8 the moderator in compartment eight, kin the turnover rate, p the amplification factor, I(Cp) the inhibitory drug function, kcap the formation of NEFA in capillaries, and kout the fractional turnover rate. For describing the inhibitory drug mechanism, the following function was used

where Imax, IC50, and γ are the maximum drug-induced inhibitory effect, plasma concentration at 50% reduction of maximal effect (potency), and sigmoidicity factor, respectively. The turnover of the moderators was given by

where ktol is a fractional turnover rate constant. NEFA build-up R and moderator M (where i = 1, …, 8) at baseline and steady state becomes M = R and the relation between turnover rate kin, the baseline NEFA concentration, R0, kout, and kcap becomes

This relationship can be used to eliminate a redundant parameter from the model, for instance by expressing kin in terms of kout and R0. The eighth moderator compartment, M8, stimulates the loss of R. The dual-action of insulin on NEFA regulation is captured by M1 and M8, where M1 describes the rapid inhibition of the hydrolysis of TG to NEFA and glycerol, and M8 the delayed stimulation of reesterification of NEFA to TG.19–22 The moderator was affected by R via a first-order process, ktol R, and each transduction step used the same transit time, 1/ktol. When NiAc inhibits the formation of R, it decreases. Consequently, the production of moderator M1 will decrease. Since the formation of NEFA is inversely proportional to the moderator raised to the power of p, the formation of NEFA will thus increase. After a delay, the level of moderator in the final compartment, M8, will also decrease, resulting in a smaller loss of NEFA. Eventually, the concentrations of R and M (i = 1, …, 8) will equilibrate.

Table 2

Initial Parameter Estimates for Pharmacodynamic Model of NEFA Plasma Concentration

Parameter	Definitions	Normal Sprague–Dawley Rats	Obese Zucker Rats
R0 (mmol L−1)	Baseline NEFA concentration	0.2	0.05
kout (L mmol−1 min−1)	Fractional turnover rate	0.03	0.03
ktol (min−1)	Turnover rate of moderator	0.001	0.001
kcap (mmol L−1 min−1)	Turnover rate of moderator	1	1
p	Amplification factor	0.08	0.08
IC50 (μmol L−1)	Potency	1	1
γ	Sigmoidicity factor	1	1

Joint Feedback Model of NEFA for Both Normal Sprague–Dawley Rats and Obese Zucker Rats

Pharmacodynamic parameters have previously been shown to differ between normal and obese rats.4,5 It is therefore important to quantify these differences as this information could be used to scale from the normal model to the disease model in, for instance, a compound ranking situation. In previous studies, the main differing characteristics were an increased NEFA baseline, R0, more distinct development of tolerance during NiAc administration and a lower rebound peak.3–5 The parameter estimate that differed the most except the baseline estimate R0 was the fractional turnover kout. However, these are functionally related with kin which is a derived parameter in our model see Eq. (10). Also, there was a noticeable difference in γ, the sigmoidicity factor.

Our suggestion is to design a covariate model to handle both normality and disease in the NiAc/NEFA model. The parameters identified to differ between the groups were replaced with a categorical variable dependent on whether the subject belonged to the obese or normal group. For example, for the NEFA baseline R0, we have the modified parameter for individual i expressed as

where Z is either 1 or 0 depending on whether the individual belongs to the normal or obese group, respectively, and R0N is the parameter value for the normal animals and R0D for the diseased animals. Equivalent expressions are used for the parameters kout and γ. For every disjoint parameter between normal Sprague–Dawley rats and the obese Zucker rats this leads to at least one additional parameter to estimate, since both the actual parameter and its variation may need to be estimated. In addition, other parameters could be set to differ between the groups. As a first attempt, we chose R0, kout, and γ to find a model that can translate between normal and obese rats.

Initial Parameter Estimates

This procedure has been described previously.3,4,23 All initial estimates are shown in Table2. The parameter kin is a derived parameter, see Eq. (10). The efficacy Imax is fixed to 1 and the sigmoidicity, γ, and amplification factor, p, are initially set to 1. The parameters R0, kout, kcap, ktol, and IC50 are all derived manually from data. In response to a high dose of NiAc (C » IC), Eq. (5) is approximated by

Provided kcap is initially much less than kout RR0, Eq. (12) can be simplified to

where M8 is approximated by R0 because of the delay caused by the cascade of moderators. Thus, the initial downswing of R on a semi-logarithmic plot gives a slope of −koutR0. The lower physiological limit of NEFA was reached following the highest dose of NiAc. Hence, kcap was estimated from

where Rss denotes the manually approximated steady state response. The initial estimate of ktol can be approximated from the time course of the log-linear decline of NEFA postrebound. The IC50 value was approximated by the concentration of NiAc resulting in a half-maximal response following a 30 min infusion.

Parameter Estimation

Parameters are estimated by optimizing the approximate population likelihood function resulting from the FOCE approximation, see the Supporting material. The optimization problem is solved using the gradient-based method BFGS.24 The gradient of the objective function needed by such methods is typically computed by finite difference approximations. However, finite difference approximations might become an unreliable description of the gradient due to the numerical solutions of the model equations. Numerical ODE solvers using adaptive step length are known to introduce quantification errors to the objective function, making it non-smooth on small scales.25,26 To overcome such problems, the gradient can be determined by formally differentiating the objective function Eq. 26. To calculate the components of the gradient w.r.t. the parameters, one typically needs to solve the so called sensitivity equations.27,28 These equations are obtained by differentiating the system equations with respect to the parameters to be estimated. Not only is the approach of using sensitivity equations more accurate, it is also generally faster since a lot of likelihood function evaluations are avoided. The estimation routine was implemented in Mathematica 9 (Wolfram Research, 2012).

Mixed-effects modeling of the NiAc/NEFA system were previously performed separately for normal rats and obese rats by utilizing NONMEM (Version VI level 2.1; Icon Development Solutions, Hanover, Maryland).3,4,23 Interindividual variability was modeled as lognormally distributed parameters for all disposition parameters of NiAc and NEFA and the random residual variability was modeled as a function of proportional or additive error for the PK and PD. The individual disposition parameters of NiAc were introduced as fixed parameters in the analysis of NEFA data.

Results

Disposition Analysis of NiAc

The parameters of the PK model have previously been estimated3 using time series data of NiAc plasma concentration. Hence, these values will be used in present analysis. The parameter values are shown in Table1.

Table3 displays the resulting parameter estimates together with the results from our previous study.4 The parameter estimates are similar to those previously found. However, our method gives a more robust estimate in terms of initiation of the estimation. Using NONMEM, it is possible to reach the same optimum by manually restarting the estimation algorithm closer and closer to this point in the parameter space and acknowledging that the approximate population likelihood function becomes larger. With our more accurate method for calculating the gradient, we can start the estimation in the derived initial values and still reach this set of parameter values. This gives a less complex and subjective optimization.

Table 3

The Ratios Between Estimated Parameters in the Joint and Separate Models

koutD/N	1.18
γD/N	0.138
R0D/N	1.443

The observed time profiles of NEFA plasma concentrations for all doses and subjects, including both obese Zucker rats and normal Sprague–Dawley rats, are shown together with the results of the population model simulations using the estimated parameters in Figure 3. In all cases, NiAc administration decreased NEFA plasma concentrations. As NiAc concentration is decreased there is a rapid return to, and above, the predose baseline concentration of NEFA. An oscillatory behavior is observed post rebound until the system has returned to baseline. The magnitude of the overshoot increases with increasing NiAc exposure, either in duration of administration or amount. Simulated responses of a typical individual (with random effect parameters set to zero) are shown together with the population variability bands calculated by the Monte Carlo method.29 The simulated variability bands give additional information about the variability in population dynamics. The variability bands were constructed in the following way. First, a large number of parameter sets were sampled according to their estimated distributions. Then, each parameter set was used for simulating the PD model. The bands mark, at each simulated time point, the 5% and the 95% quantile of the simulated NEFA plasma concentrations. Note that complete bands do not correspond to a model simulation with a specific parameter set, but are a representation of all the 90% intervals at each simulated time point. Comparing the experimental data to the variability bands revealed that 1274 data points out of in total 1440, or 88.5%, resided within the bands.

Figure 3

Observed NEFA plasma concentration–time profiles together with population model fits of NEFA plasma concentrations–time data for normal Sprague–Dawley rats after infusion (a–h), oral administration (i–k), for obese Zucker rats after infusion (l–m). Infusion of (a) 1 μmol kg−1; (b) 5 μmol kg−1; (c) 20 μmol kg−1, over 30 min; (d) 5 μmol kg−1; (e) 10 μmol kg−1; (f) 51 μmol kg−1, over 300 min; (g) 5 μmol kg−1 over 30 min followed by a stepwise decrease in infusion rate every 10 min for 180 min; (h) 5 μmol kg−1 over 30 min followed by a stepwise decrease in infusion rate every 10 min for 180 min, and another 5 μmol kg−1 infusion over 30 min. Oral dose of (i) 24 μmol kg1, (j) 82.1 μmol kg−1, and (k) 812 μmol kg−1. Infusion of (l) 20 μmol kg−1 over 30 min and (m) 51 μmol kg−1 over 300 min. Shaded bands show the Monte Carlo 90% population variability bands that describe the variability of the population.

The model also performs well at the individual level, showing high consistency between experimental data and model predictions. This is illustrated in Figure 4, where representative individual fits are shown for all administrations. Because of its flexibility the model captures, the vast differences of the various drug provocations. However, the tolerance build up for the longer infusion experiments is not completely captured, which might indicate a possible need for an additional mechanistic component in the model.

Figure 4

Representative individual model fits together with individual data of NEFA plasma concentration–time data for normal Sprague–Dawley rats after infusion (a–h), oral administration (i–k), for obese Zucker rats after infusion (l–m). Infusion of (a) 1 μmol kg−1; (b) 5 μmol kg−1; (c) 20 μmol kg−1, over 30 min; (d) 5 μmol kg−1; (e) 10 μmol kg−1; (f) 51 μmol kg−1, over 300 min; (g) 5 μmol kg−1 over 30 min followed by a stepwise decrease in infusion rate every 10 min for 180 min; (h) 5 μmol kg−1 over 30 min followed by a stepwise decrease in infusion rate every 10 min for 180 min; and another 5 μmol kg−1 infusion over 30 min. Oral dose of (i) 24 μmol kg−1, (j) 82.1 μmol kg−1, and (k) 812 μmol kg−1. Infusion of (l) 20 μmol kg−1 over 30 min and (m) 51 μmol kg−1 over 300 min.

Figure 5a shows the concentration–response relationship at steady state for the separate models. It shows an upward shift and a more shallow behavior in the concentration–response relation for the obese animals compared to the normal ones. The normal profile is also more s-shaped than the disease profile which is flatter. Population variability bands, resulting from the interindividual variability of the parameters, are displayed. Note that the bands do not correspond to a response–concentration relationship for an actual individual but to the population variability in responses at each simulated concentration.

Figure 5

Simulated steady-state plasma concentration versus predicted NEFA plasma concentration at equilibrium with population variability bands. These results were derived using the estimated parameters from Table3 for obese (solid) and normal rats (dashed) estimated separately (left) and from the joint (right) model at steady state. Red bands correspond to obese rats and blue bands to normal rats.

Joint Analysis and Disease Model

The joint NEFA feedback model is defined in section Joint Feedback Model of NEFA for Both Normal Sprague–Dawley Rats and Obese Zucker Rats. The parameter estimates for this model are shown in Table3. For normal rats, kout was decreased in the joint analysis compared to the separate analysis. Conversely, the value of the parameter kout for obese rats was increased, becoming more similar to the value of the normal rats. The estimated values of the sigmoidicity factor, γ, and NEFA baselines, R0, did not differ significantly for either normal or obese rats compared with the separate analysis. Table3 also shows the values of the interindividual variation (IIV), which generally increase in the joint analysis compared to the separate analysis.

Figure 5b shows the concentration–response relationship at steady state for the joint model. The response is similar to the separate models but the population variability is increased, especially for obese rats. We also note a slightly decreased slope in the response for obese rats in the joint model.

Figure 6 shows population model fits from the separate analysis together with results from the joint analysis for a 30 min infusion of 20 μmol kg−1. Figures 6a and 6b show the results for normal animals in the separate and joint model, respectively. Figures 6c and 6d show the corresponding results for obese animals. We observe that the dynamics for the normal group does not change extensively. However, for the obese group the population model curve slightly changes shape when the information from the normal group is included. Additionally, the variability in obese response in the joint model is markedly increased compared with the separate analysis. In Table4, the ratios between the joint and separate parameter estimates are listed.

Figure 6

Population model fits of NEFA plasma concentration–time data for normal Sprague–Dawley rats after infusion (a and b) with 20 μmol kg−1 for 30 min and (c and d) the corresponding plots for obese Zucker rats. Panels (a and c) correspond to NEFA plasma concentrations simulated with the parameter estimates resulting from the separate analysis and panels (b and d) to the estimates from the joint analysis, respectively. Shaded bands show the Monte Carlo 90% population variability bands that describe the variability of the population.

Table 4

Population Pharmacodynamic (NEFA Response) Parameter Estimates and Interindividual Variability (IIV) with Corresponding Relative Standard Errors (RSE%)

		Joint Analysis		Separate Analysis		Ahlström et al.
Parameter	Definitions	Estimate	IIV	Estimate	IIV	Estimate	IIV
Normal Sprague–Dawley Rats
R0 (mmol L−1)	Baseline NEFA conc.	0.741(4.42)	19.7(31.3)	0.739(3.50)	19.9(28.1)	0.736(4.33)	21.6(27.3)
kout (L mmol−1 min−1)	Fractional turnover rate	0.290(10.1)	71.4(11.4)	0.316(9.83)	47.4(10.7)	0.273(10.2)	42.7(11.3)
ktol (min−1)	Turnover rate moderator	0.0245(2.00)	–	0.0239(1.92)	–	0.0231(1.90)	–
kcap (mmol L−1 min−1)	NEFA form. in plasma	0.0245(10.2)	–	0.0281(9.85)	–	0.0230(10.1)	–
kin (mmol L−1 min−1)a	Turnover rate NEFA	0.0940	–	0.102	–	0.0844	–
p	Amplification factor	1.20(4.02)	–	1.15(2.61)	–	1.13(2.76)	–
IC50 (μmol L−1)	Potency	0.0820(15.5)	129(40.1)	0.0757(14.9)	134(36.2)	0.0680(15.4)	131(34.9)
γ	Sigmoidicity factor	2.16(4.87)	–	2.19(4.33)	–	2.18(4.48)	–
Imax	Efficacy	1	–	1	–	1	–
σ1	Residual prop. error	–	–	–	–	–	–
σ2	Residual add. error	0.0110(2.84)	–	0.0102(2.45)	–	0.00913(2.63)	–
Obese Zucker rats
R0 (mmol L−1)	Baseline NEFA conc.	1.07(6.34)	24.0(40.6)	1.04(6.31)	15.8(34.2)	1.06(6.52)	14.8(34.6)
kout (L mmol−1 min−1)	Fractional turnover rate	0.343(31.1)	151(72.1)	0.0683(24.0)	29.7(63.6)	0.0986(24.5)	69.5(62.9)
ktol (min−1)	Turnover rate moderator	0.0245(2.00)	–	0.0697(9.87)	–	0.0297(10.2)	–
kcap (mmol L−1 min−1)	NEFA form. in plasma	0.0245(10.2)	–	4.94 × 10−8(−)	–	0(fixed)	–
kin (mmol L−1 min−1)	Turnover rate NEFA	0.399	–	0.0837	–	0.125	–
p	Amplification factor	1.20(4.02)	–	3.19(15.3)	–	2.01(15.8)	–
IC50 (μmol L−1)	Potency	0.0820(15.5)	129(40.1)	0.0433(49.6)	74.7(39.6)	0.0538(56.5)	135(43.4)
γ	Sigmoidicity factor	0.298(9.63)	–	0.563(8.99)	–	0.347(9.34)	–
Imax	Efficacy	1	–	1	–	1	–
σ1	Residual prop. error	–	–	–	–	12.0(19.9)	–
σ2	Residual add. error	0.0110(2.84)	–	0.0111(2.56)	–	–	–

Calculated as a secondary parameter.

The joint model uses categorical covariates to make a distinction between healthy and diseased subjects. However, we may hypothesize the existence of disease conditions lying in-between these two groups. We suggest that our joint model can be used to describe various “degrees of disease” by letting the variables Zi take values in the continuous range between 0 and 1. Based on this idea, we decided to quantify the extent and variability of the rebound effect following the termination of a 300 min infusion. The extent of rebound was defined as the maximum level of NEFA reached, expressed in percentage above the baseline. Simulations were performed for six different levels of degree of disease, and for three different doses, see Figure 7. Monte Carlo simulations of 1000 sampled individuals were used for each condition. Since PK models differed for normal and obese rats, both PK models were solved simultaneously and the solutions were weighted accordingly. For the two lowest doses, the median extent of the rebound was highest for intermediate degrees of disease, whereas the highest dose showed a monotonically decreasing median rebound with increasing degree of disease. For all doses, there is a clear trend of markedly decreased variability in the extent of rebound as function of increased degrees of disease.

Figure 7

The effect of degree of disease on extent and variability of rebound. The extent of rebound following the termination of a 300 min infusion was quantified as the maximum level of NEFA reached, expressed in percentage above the baseline. The population variability of this quantity was calculated for different degrees of disease (0%, 20%, 40%, 60%, 80%, and 100%, respectively) for three different doses (5, 10, and 51 μmol kg1, respectively). Results are presented as box and whisker plots, showing the median values, 25% and 75% quantiles, respectively, and the minimum and maximum observations. For each dose, a trajectory is joining the medians to emphasize the changes in extent of rebound as a function of the degree of disease.

Sensitivity Analysis

Sensitivity analysis is used to determine which parameters are most important in a model.30 It can answer questions about which parameters that affect a measurable entity the most, or how easy it is to identify one specific parameter compared with others. This information may be used for experimental design, including when and how to take samples. Sensitivity analysis is performed by analyzing how a small variation of a parameter value relates to a small variation of a response variable.

We performed this analysis for all the parameters of the PD system with the NEFA plasma concentration as the response variable. Figure 8 shows a panel plot of the time-dependent sensitivities for one NiAc administration route, namely constant infusion for 300 min with a dose of 51 μmol kg−1 of NiAc. The plots are ordered in size of magnitude for the derivatives of the response w.r.t. each parameter. The parameter ktol affects the response the most of all the PD parameters, followed by kcap.

Figure 8

Time-dependent sensitivity analysis for a constant infusion for 300 min with a dose of 51 μmol kg−1 of NiAc in normal rats. The response, plasma concentration of NEFA, has been differentiated with respect to each parameter of the pharmacodynamic system.

Discussion

Compound ranking in early drug discovery phase is often based on limited results obtained from animal models of normal or special populations. It is therefore important to understand and quantify the PD differences between disease and normality. Obesity causes physiological changes that can alter how a drug regulates turnover of the response.6–10 In particular, NEFA turnover is altered to respond differently to NiAc in the obese Zucker rat compared with the normal Sprague–Dawley rat. Therefore, the amount of NiAc required to reach the desired effect is expected to be different in normal and obese animals.

Previous studies of the differences between normal Sprague–Dawley rats and Zucker rats that were not obese showed that they have similar lipid and lipoprotein profiles and small differences in physiological characteristics such as insulin levels.31,32 In contrast are the obese Zucker rats, where the physiological characteristics differed significantly to the lean Sprague–Dawley rats, with higher insulin concentrations and NEFA baseline was approximately 70% higher than for the lean ones. The obese rats showed insulin resistance since blood glucose levels were the same as in normal rats in spite of the fact that insulin concentrations were significantly higher in obese animals. In this study, we further evaluated and quantified the difference in PD between normal Sprague–Dawley rats and obese Zucker rats.

Administration of NiAc led to a rapid decrease in NEFA plasma concentrations, and previously published results show that these concentrations have a lower physiological limit.3–5 During longer infusions, we observe tolerance, by which the NEFA concentrations slowly rise with a rapid rebound above the predose baseline when the infusion is terminated. Tolerance build-up occurs in response to all types of NiAc administrations, but was less evident for obese Zucker rats. These rats have developed a resistance to insulin and their NEFA response to changes in NiAc concentration is less sensitive than that of the normal rats. The tolerance builds up more smoothly and the overshoot is not as pronounced as for the normal rats. The discrepancy between normal and obese rats leads us to the assumption that disease changes physiological effects on the PD. Thus, the joint model can be an appropriate tool for translating from the normal rat model to the disease rat model when such translation is necessary.

Our investigation of the NEFA PD showed that, the main differences between obese and normal rats were an increased NEFA baseline concentration, more distinct development of tolerance during NiAc administration and a lower rebound peak. Parameter estimates that differed the most were the baseline estimate R0, the turnover rate kin, and the fractional turnover kout. These are functionally related, as indicated in Eq. (10). The more flat and upwardly shifted relationship between NiAc concentration and NEFA response explains the greater extent of tolerance and lack of rebound in obese animals. Also there is a pronounced difference in the sigmoidicity parameter, γ. In conclusion, these results are similar to previous results.4

The population variability bands introduced in our study account for the variation in system dynamics that is due to interindividual parameter variability. This is an important extension to previous results, which did not address variations in dynamics and response. However, our analysis is restricted to the assumption of uncorrelated random parameters. This assumption gives a tractable optimization problem in terms of parameter identifiability, but may lead to an overestimation of the population variability.

Concentration–Response Relationship

Figure 5a shows the concentration–response relationship between NiAc and NEFA at steady state for normal (dashed) and obese (solid) for the separate models. Comparing the normal and obese responses shows, as expected, an upward shift and a more shallow concentration–response relation for the obese animals. The normal profile is also more s-shaped than the disease profile which is much flatter. This profile is conclusive in a system with (more) insulin insensitivity, as is the case for the obese Zucker rats.4 The concentration–response curves provide a prediction of drug effects when scaling from normal individuals to individuals with disease. For instance, increasing the NiAc concentration beyond 1 μmol/L does not change the NEFA steady-state response in normal rats. However, for obese rats the response continues to decrease as NiAc steady-state concentration is further increased for another two orders of magnitude.

The joint model, with three parameters differing between the normal and the diseased group manages to capture the features separate models, as shown in Figure 6. These three parameters do give a compact description of the main characteristics that differ between normal and obese rats. Firstly, inspecting the time series data in Figure 3, it is obvious that the baseline NEFA plasma concentration, determined by the parameter R0, differs significantly between the groups. Secondly, the nonlinearity of the concentration–response relationships is a consequence of the sigmoidicity parameter γ. The dramatically smaller γ in diseased animals suggests an inherent resistance of the system to changes in NiAc concentration. Finally, since the down-swing of NEFA differs between the two groups, one also expects the fractional turnover rate kout to differ. Comparing the parameter values between the two groups seen in Table4 as ratios gives an idea of how to quantify the magnitude of disease by concentrating the information from the parameter estimation into a more compact description. The results from the joint-analysis build on the strength of the extensive data from normal rats. The uncertainty of the IC50 value for obese rats is reduced fourfold and we gain significant higher precision in our parameter estimates. The combination of fast NiAc kinetics with a build-up of tolerance in the NEFA response can lead to extensive rebound effects. This is undesirable since NEFA levels much higher than the baseline level is a potential safety issue. It has previously been shown that the extent of rebound is influenced by the duration of NiAc exposure and that rebound is substantially higher in normal rats. Using the joint model derived in this study, we have further characterized the rebound effect by simulating its population variability and investigating its dependence on a hypothesized degree of disease (Fig. 7). Interestingly, for lower doses it appears that the rebound effect may be strongest for intermediate degrees of disease. We additionally found that the population variability of the rebound effect has a strong dependence on the degree of disease, with comparatively little variability among the fully diseased animals. Taken together, these simulation results suggest that the rebound effect, and hence its impact on the dose planning for an animal disease model, may be very sensitive to the exact degree of disease.

In conclusion, joint models can be used as a tool when screening drugs on normal animals and then translating the results to the disease model, giving scientists a way to quantify and evaluate the differences between the normal subject and the diseased subject.

New Estimation Algorithm

The parameters of the present analysis were estimated by optimizing the approximate population likelihood function of the problem resulting from the FOCE method but with gradients computed by a method that has higher precision than the finite difference approximation of the gradient used in NONMEM.13 Because of the higher precision, we were able to provide better estimates for the IIV of IC and kout, which was shown to be significantly decreased compared with previous work. We also note that our method, based on the sensitivity equations, can reach the optimum from the manually derived initial estimates without repeatedly restarting the optimization from new initial values. Compared to NONMEM, this gives a more robust calculation of the optimal parameter set in the sense of convergence properties since no manual supervision is needed.

Some of the sensitivity profiles in Figure 8 showed similar behavior, for example, IC50 and kin and, γ and p. These similarities indicate dependencies between parameters, which can lead to identifiability issues in the parameter estimation. For future models, it may be valuable to reduce or reparameterize the model in order to improve identifiability. Sensitivity analysis also shows which parameters are important at different transient behaviors. We deduced that many of the parameters affect the response during the onset and the end of NiAc administration, and during the rebound phase. For instance, the fraction turnover rate ktol is the parameter that affects the NEFA plasma concentration response the most. Similarly, kin, the turnover rate of NEFA, is most important at the end of infusion.

The analysis is carried out separately for each parameter keeping all but one parameter fixed. The ultimate goal would be to also include expected or estimated covariation of parameters in the sensitivity analysis to better characterize the actual system response to observed variation.

Another use of sensitivity analysis is to guide the design of experiments. The sensitivity profiles in Figure 8 show the importance of dense sampling at the start and end of infusion, when the response is most rapid. It is common to have fast transients in the beginning of an administration, and in some studies it is known to be valuable to sample often at the start. However, our analysis shows that the parameter estimates are substantially affected by information collected post infusion. We suggest a design, if possible for ethical reasons, with more samples taken at the end of the experiment for increased precision of the parameter estimates.

Recommendations and Value for Drug Discovery

Disease may alter both the PK and the PD of a drug, resulting in an altered concentration–response relationship, resulting in an altered concentration–response relationship. Therefore, a quantitative approach in modeling is important even in the early stages of drug discovery.

Administration of NiAc leads to a tolerance build-up and to extensive rebound in normal Sprague Dawley rats, but not in obese Zucker rats.4,5 The high rebound can become a safety issue. Therefore, the population variability bands that we obtained provide important knowledge about the variation in individual response.

Conclusions

In conclusion, our study shows that both the separate and the joint model capture the PD characteristics successfully. The joint model describes how normal and obese rats differs with only three separate parameters, thereby yielding important information for predicting outcomes in diseased subjects based on data from normal ones. Another important conclusion is that the planning of experiments and dose scheduling cannot rely on pilot studies on normal animals alone. The obese rats have a much less pronounced concentration–response relationship and need much higher doses to exhibit desired response. Building joint normal-disease models with scaling parameter(s) to characterize the “degree of disease” can be a very useful tool when designing informative experiments on disease animals.