Investigating Stochastic Differential Equations Modelling for Levodopa Infusion in Patients with Parkinson's Disease.

BACKGROUND AND OBJECTIVES: Levodopa concentration in patients with Parkinson's disease is frequently modelled with ordinary differential equations (ODEs). Here, we investigate a pharmacokinetic model of plasma levodopa concentration in patients with Parkinson's disease by introducing stochasticity to separate the intra-individual variability into measurement and system noise, and to account for auto-correlated errors. We also investigate whether the induced stochasticity provides a better fit than the ODE approach.
METHODS: In this study, a system noise variable is added to the pharmacokinetic model for duodenal levodopa/carbidopa gel (LCIG) infusion described by three ODEs through a standard Wiener process, leading to a stochastic differential equations (SDE) model. The R package population stochastic modelling (PSM) was used for model fitting with data from previous studies for modelling plasma levodopa concentration and parameter estimation. First, the diffusion scale parameter (σw), measurement noise variance, and bioavailability are estimated with the SDE model. Second, σw is fixed to certain values from 0 to 1 and bioavailability is estimated. Cross-validation was performed to compare the average root mean square errors (RMSE) of predicted plasma levodopa concentration.
RESULTS: Both the ODE and the SDE models estimated bioavailability to be approximately 75%. The SDE model converged at different values of σw that were significantly different from zero. The average RMSE for the ODE model was 0.313, and the lowest average RMSE for the SDE model was 0.297 when σw was fixed to 0.9, and these two values are significantly different.
CONCLUSIONS: The SDE model provided a better fit for LCIG plasma levodopa concentration by approximately 5.5% in terms of mean percentage change of RMSE.

Key Points

Introduction

Parkinson’s disease is the second most common neurodegenerative disease after Alzheimer’s disease [1]. While there is no cure for Parkinson’s disease, levodopa medication is widely used as the ‘gold standard’ for treatment [2]. Treatment in advanced stages can be provided by oral tablets or intestinal administration of levodopa/carbidopa gel (LCIG) through an adjustable pump. LCIG administration reduces the variation in levodopa concentration [3, 4].

Chan et al. [5] suggested a two-compartment pharmacokinetic model for levodopa concentration over a one-compartment model, as the two-compartment model best described their intravenous data. Westin et al. [6] further developed this model to include the pharmacodynamic part for duodenal LCIG infusion and described it by four ordinary differential equations (ODEs) and a sigmoid effect equation. Simon et al. [7] developed a one-compartment combined pharmacokinetic/pharmacodynamic model in patients receiving levodopa orally. Othman and Dutta [8] also found that the two-compartment model fits data better for both oral and infusion administration.

In this study, the two-compartment model for infusion developed by Westin et al. [6] is considered, where the first three equations comprise the pharmacokinetic part, which describes the plasma levodopa concentration of LCIG in the system after dose and/or infusion. In the study by Westin et al. [6], the absorption parameters were estimated from dose LCIG plasma levodopa concentration, and the pharmacodynamic parameters from concentration–time profile and effect data.

The model developed by Westin et al. [6] does not include covariates other than weight that could affect LCIG levodopa pharmacokinetics, so they are included in the inter-individual and intra-individual variability. A dynamic system such as this is driven both by our own control inputs and disturbances which is out of our control and cannot be modelled deterministically [9]. Furthermore, measurement tools introduce their own system dynamics and distortions. These factors necessitate the separation of the sources of variability.

Furthermore, pharmacokinetic/pharmacodynamic models such as the one developed by Westin et al. [6] have generally been based on ordinary differential equations (ODEs) by making assumptions of independence or non-correlation between the residuals. Correlation between residuals is not that uncommon, however, and may lead to inaccurate estimates of the parameters, as demonstrated by Karlsson et al. [10]. According to Overgaard et al. [11], models accounting for correlated residuals also produce better estimates of the inter-individual variation and the structural parameters.

The intra-individual variations may be more precisely modelled and the residual correlation could be accounted for by using stochastic differential equations (SDEs) [11, 12]. SDEs can be used to explain the discrepancies between individual predictions and observations by separating the noise into dynamic system noise and measurement noise [11]. The dynamic system noise comes from the dynamics of the system itself and could arise from model deficiencies or true random fluctuations in the system. The measurement noise, on the other hand, comes from the uncorrelated part of the residual variability. The system noise allows the modeller to compensate for unknown phenomena such as modelling errors, approximations, and oversimplifications [11, 13]. However, there is a complication when fitting SDE models because optimisation algorithms often fail to find the global optima, because “the stochastic volatility term roughens the likelihood surface and thereby inflates the fit of many different parameter sets that would otherwise not have adequately accounted for the data” [14].

Stochastic pharmacokinetic/pharmacodynamic models have been demonstrated in other areas. Donnet and Samson [15] reviewed estimation methods of SDEs for various pharmacokinetic/pharmacodynamic models. Leander et al. [13] performed a simulation study of a stochastic pharmacokinetic model and studied the application on nicotinic acid disposition in obese Zucker rats with a one-compartment model. Klim et al. [16] developed the R package population stochastic modelling (PSM) for mixed-effects models based on SDEs, and demonstrated an example with insulin secretion rates on a two-compartment model.

It is important to individualise the pharmacokinetic profile of patients and obtain precise predictions of the distribution of levodopa in their body, which in turn allows more refined predictions of the effect. Therefore, the purpose of this study is to investigate a two-compartment pharmacokinetic model of levodopa in patients with Parkinson’s disease by introducing stochasticity to account for correlated residual errors and check whether the SDE model fits better with the data than its ODE counterpart through cross-validation.

Theory

In this section the stochastic mixed-effects theory is briefly introduced. First the ODE model is presented followed by the SDE model, where a stochastic Wiener process is added (in the mathematical expressions, boldface letters indicate vectors throughout).

ODE Model

Non-linear mixed-effects models, with ODE, are used to describe data with the following general structure:where y is a vector of measurements or observations at time t for individual i, N is the number of individuals and n is the number of measurements for individual i. In such a model the variation is separated into intra-individual and inter-individual variations. For more details please refer to Refs. [11, 17].

The intra-individual variation is given by:where t is the continuous time variable, xi(t) is the state of the model at time t (such as amount or concentration of drug), ui(t) is a vector of optional inputs (such as dose administered or infusion), t is the jth time measurement, and e is the measurement error, where e ~ N (0, ). The distributions of the residuals for each subject were visually inspected and found to be normally distributed for the majority of the subjects, so an additive error model (3) was chosen. Furthermore, i is a vector of parameters for each individual i including both population level fixed effects and individual-specific random effects i. Hence, i captures the inter-individual variation.

SDE Model

The ODE model may be extended to an SDE by including an additional source of variation in the first-stage model, the system noise. This separates the intra-individual variation into two sources of noise—(1) measurement noise and (2) system noise, driven by a newly introduced parameter matrix σ. The system noise, which is a continuous stochastic process, is added to the differential equations to allow for random variations in the evolution of the states. The measurement noise, on the other hand, enters the model at discrete time intervals and is included in the measurement equation [13]. The extended first-stage model then becomes:where σwdw is the stochastic part of the model and is referred to as the system noise. It should be noted that if the magnitude of the system noise σw is assumed to be zero, then the SDE (4) simply reduces to the ODE (2). The component dw refers to the infinitesimal increments in the noise process (). This individual system noise is a standard vector Wiener process, which is a continuous time Gaussian process such that the mean and variance between any two time points are:where I is the identity matrix. The dw Wiener increments are assumed to be independent across individuals and independent of the measurement error.

The SDE approach offers more flexibility by introducing an error in the differential equations, as opposed to the ODE framework where the error is only in the measurement equation. It may also account for auto-correlated residuals and allow for modelling different residual error correlation patterns [11].

With the updated model, variability is now separated into three sources—measurement noise, system noise, and inter-individual variability. To fit this model, the R package PSM [16] was applied, which uses a Kalman filter to solve the system of SDEs. For more information, please refer to the Mathematics Guide to CTSM [17] or the original work by Kalman [18].

Methods

Studies and Subjects

The data used in this study were pooled from two of the three studies investigated by Westin et al. [6] (referred to as study 1 and study 2 throughout). Data from their first study showed LCIG plasma levodopa concentrations that were generally stable, so were not used for our study. Some of the model parameters had been fixed to values found in the literature [5]. In our study, the parameters are fixed according to the aforementioned literature [5, 6]. The baseline characteristics for the patients in the studies used are given in Table 1.

Table 1

Baseline characteristics for the patients in the reference study by Chan et al. [5] and two other studies [6, 19] used for modelling (mean ± standard deviation)

Characteristic	Chan et al. [5]	Study 1 [6]	Study 2 [19]
No. patients (male/female)	20 (12:8)	3 (3:0)	5 (3:2)
Weight, kg	78.7 ± 12.4	69.7 ± 14.2	62 ± 8.4
Age, year	59.8 ± 10.7	62.3 ± 2.5	60.8 ± 6.1
Age at onset, year	56.1 ± 10.9	44.3 ± 5.7	48.8 ± 7.9
Duration of disease, year	3.8 ± 1.7	18.0 ± 4.6	12.0 ± 2.9
Hoehn and Yahr stage	2.5 ± 0.6	2.7 ± 0.6b	4.4 ± 0.9a

aStage assessed at worst

bStage assessed during ‘on’

Study 1 included three patients who had advanced Parkinson’s disease, were already on Duodopa®, and had no concomitant diseases (the data were from 2 nonconsecutive days for each patient) [6]. As one of the patients did not attend the second day, the data are from five different occasions. During Study 1, the patients were given a bolus dose in the morning with the LCIG. When the patient’s clinical effect had reached baseline again, the patient’s normal infusion rate was started and data were collected for 2 h. Blood samples and effect measurements were collected every 5–10 min. The plasma samples were analysed by high-performance chromatography as described by Nyholm et al. [3]. No interacting drugs were administered to the patients during the study. This study protocol was accepted by the ethics committee of the Karolinska Institute, Sweden, and the patients had given informed consent in accordance with the Helsinki declaration.

Study 2 included five patients with idiopathic Parkinson’s disease who were difficult to keep in ‘on’ state without dyskinesia (the data were from three occasions for each patient) [19]. Data were collected from the patients for five 4 h-periods with five different infusion rates over 2.5 days. On day 1, doses that were used in the first half-day were well-adjusted and optimised for each individual. For day 2, 120% of the optimised dose was used for the first half-day and 90% of the optimised dose on the second half-day. For day 3, 80% of the optimised dose was used for the first half-day, and 110% of the optimised dose on the second half-day. The doses were blinded on days 2 and 3. Data were collected every 20–30 min. The total number of occasions for study 2 was 15, and no interacting drugs were administered to the patients during the study. The study was approved by the local ethics committee and the Swedish Medical Products Agency.

Modelling

The pharmacokinetic model of levodopa is a system of three differential equations (Eqs. 8–10) as described by Westin et al. [6]. Initially, the ODE-based pharmacokinetic model was fitted with all 20 occasions and the residuals were visually checked for correlation for each patient.

It is assumed that three sources of variability (measurement noise, system noise, and inter-individual variability) may be separated. A system noise σwdw is introduced in the second equation (Eq. 9), which is for the second compartment and where the observation is made, to extend the model into stochastic partial differential equations. This choice is made because observation is taken from the second compartment, where we expect a disturbance to be created in the system and thus introduce system noise. In this model, the system noise accounts for the uncertainty in the dynamics that may arise from true fluctuations of the system, or any oversimplification of the model:where Inf represents levododa infusion (mg/min), ai is the amount (mg) in compartment i, Vi is apparent volume (L) in compartment i, Q is inter-compartmental clearance, BIO is the bioavailability (fraction absorbed), and σwdw represents system noise. The absorption time constant (min) pharmacokinetic parameter is defined as TABS = 1/. The structural pharmacokinetic model is shown in Fig. 1.

Fig. 1

Structural pharmacokinetic model adopted from Westin et al [6]. Inf levodopa infusion (mg/min); a amount (mg) in compartment i; V apparent volume (L) in compartment i; Q inter-compartmental clearance (L/min); CL clearance (L/min)

The variability in the parameters is modelled by Eqs. (11–15), where Pindividual represents a random variable with a normal, zero-mean distribution for parameter P. The covariate used here is weight (WT) in kg and is included in Eqs. 11, 13, and 14 according to theory-based allometry models [20]:

Individual random effects have been added as parameters to be estimated in PSM, so that each individual may have individualised parameters. These random effects are estimated by PSM using the PSM.smooth function. Estimation of the parameters was based on the dose provided to the patients and their measured LCIG plasma levodopa concentration data. The default optimiser, optim, which is included in the PSM package was used.

The model is specified by fixing most of the parameters to the population mean values and variance estimated by Westin et al. [6], while allowing for individual variation (random effect), which is estimated. The parameter values that have been fixed are shown in Table 2.

Table 2

Parameters fixed from Chan et al. [5] and Westin et al. [6]

Parameter	Population mean	Variability
TABS (min)	28.5	0.42
BIO	Estimated	Estimated
V1 (L)	11	0.1936
V2 (L)	27	0.0625
CL (L/min)	0.52	0.0729
Q (L/min)	0.58	0.2304

TABS absorption time constant, BIO bioavailability, V apparent volume in compartment i, CL clearance, Q inter-compartmental clearance

The bioavailability population parameter along with its variability are estimated as a first step by specifying the σw component to be zero (ODE model), fixing measurement noise variance to 1, and using all 20 occasions. Bioavailability is the fraction of levodopa that enters circulation after being introduced into the body.

In the second step, the σw component and measurement noise variance are estimated along with the bioavailability population parameter by using all 20 occasions, while fixing the variability of bioavailability to 0, as it was found to be too small in the first step. It is checked whether the stochastic component is significantly different from 0.

Next, to investigate whether the model with the non-zero stochastic component provides a better fit, cross-validation is performed on both ODE and SDE models. The data were divided into training sets and test sets for each patient. For the ODE model without stochasticity, the σw component is set to zero and bioavailability is estimated. For the SDE model, the σw component is set to values from 0.1 to 1 (1 being the upper bound as measurement noise variance is fixed to 1), and the bioavailability parameter is estimated. The random effects are estimated and extracted, which are fixed in the test set to individualise the parameters.

A flow chart of the cross-validation process is shown in Fig. 2. For study 1, as the patients had two occasions only, step 3 (Fig. 2) was not needed for them, and the steps were not performed on the patient with only one occasion. For the test sets, using the dose information for that occasion and the individual parameters obtained from the training set, smoothed values are calculated by using the PSM.smooth function; to obtain smoothed estimates at time t, measurements from both before and after time t are used by the PSM.smooth function [16].

Fig. 2

Steps in the cross-validation process. Study 1 has two patients with 2 occasions each, so Step 3 is not needed for those patients; Study 1 has one patient with only one occasion, so the steps are not performed for that patient. Study 2 has five patients with 3 occasions each, so all the steps are performed for each patient in this study. σ diffusion scale parameter, RMSE root mean square error

The root mean square error (RMSE) is calculated for each test set and the average of all the RMSEs was calculated for both ODE and SDE models. The average estimate and the standard error (SE) of the bioavailability were also calculated at each σw value.

To check whether there was a significant difference between the highest and lowest average RMSEs, a binomial test was performed. The difference in RMSE in each trial (at the two σw values where the highest and lowest average RMSEs are observed) was calculated and the number of positive differences are counted as successes. The total number of trials in each cross-validation routine is 34. The null hypothesis is that the probability of a success is 0.50, and the alternate hypothesis is that the probability of a success is not 0.50.

The percentage difference of the RMSE at each trial (at the two σw values where the highest and lowest average RMSEs are observed) was calculated, and the mean percentage difference was found to give a measure of difference in the predictive ability of the two models.

The R scripts used during the current study are available from the corresponding author on reasonable request.

Results

The residuals from the ODE model for each subject over time are shown in Fig. 3. Most of the subjects show some correlation in their residuals.

Fig. 3

Residual plots over time along with best fit line for ODE model of LCIG plasma levodopa concentrations for the 3 patients in study 1 (panels 1011–1031) and the 5 patients in study 2 (panels 2011–2053). In the panel number, the second last digit represents the subject and the last digit the occasion (day)

The bioavailability parameters in the ODE- and SDE-based models, and the σw parameter and measurement noise variance in the SDE model are estimated by the PSM.estimate function by using all 20 occasions in the dataset. The maximum likelihood estimates of the parameters and their 95% confidence intervals are shown in Table 3.

Table 3

Maximum likelihood estimate (MLE), standard error (SE), and 95% confidence interval (CI) using all 20 occasions

Parameter	MLE	SE	CI
BIO (ODE model)	0.756	0.043	(0.671, 0.841)
Variability of BIO (ODE model)	0.001	0.000	(− 0.001, 0.002)
BIO (SDE model)	0.753	0.040	(0.674, 0.833)
σw (SDE model)	0.099	0.009	(0.081, 0.118)
Measurement noise variance (SDE model)	1.000	0.130	(0.746, 1.254)

BIO bioavailability, ODE ordinary differential equations, SDE stochastic differential equations, σ diffusion scale parameter

It is seen that the bioavailability population parameter was estimated to be approximately 75.6% with an SE of 0.043 in the ODE model, and about 75.3% with an SE of 0.040 in the SDE model, which means that, on average, 75.6% and 75.3% of the dose will be absorbed, based on the ODE and SDE model, respectively. The variability of bioavailability was found to be 0.001 with a 95% confidence interval of (− 0.001, 0.002); as it is too small, the variability of bioavailability was fixed to 0 for the subsequent analyses. The σw component was estimated to be approximately 0.099 with a 95% confidence interval of (0.081, 0.118), which is significantly different from 0 at 5% level of significance, implying that the variance in the model has been split with a Wiener noise component. The measurement noise variance was estimated to be 1.000 with a 95% confidence interval of (0.746, 1.254). Thus, for the cross-validation routines, the measurement noise variance was fixed to 1 to reduce computation time.

As a visual check of how well the two models fit, the smoothed estimates are plotted against the actual measurements in the same plot for each occasion. As the model fits are not distinguishable visually, they are not presented here.

Table 4 shows the average of RMSEs, bioavailability estimates, and SE of bioavailability for each cross-validation. A plot of the average RMSE from the cross-validations against the respective σw values is shown in Fig. 4. It can be seen in Table 4 that the mean of the bioavailability estimates from all cross-validations at each σw value is generally estimated to be around 75%. The average SE of the bioavailability estimates seems to have an increasing trend with increasing σw values. The largest average RMSE is 0.313 when σw is fixed to 0 (ODE model). The lowest average RMSE is 0.297 when σw is fixed to 0.9 (SDE model).

Table 4

Averages of root mean square error (RMSE), bioavailability (BIO), and standard error (SE) of bioavailability from cross-validations at different diffusion scale parameter (σw) values

σw	Average RMSE	Average BIO	Average SE of BIO
0.0	0.313	0.757	0.043
0.1	0.311	0.752	0.037
0.2	0.310	0.758	0.046
0.3	0.308	0.748	0.036
0.4	0.309	0.759	0.038
0.5	0.308	0.758	0.042
0.6	0.303	0.749	0.037
0.7	0.299	0.747	0.045
0.8	0.298	0.757	0.049
0.9	0.297	0.759	0.042
1.0	0.307	0.761	0.048

Fig. 4

Average root mean square error (RMSE) from cross-validation study for each diffusion scale parameter (σw) value

Comparing the RMSEs at each trial for these two σw values, it was found that 33 of the 34 trials yielded a lower RMSE for the SDE model compared to the ODE model. The p value for the binomial test is almost 0, therefore at the 95% confidence level, these two average RMSE values are significantly different. The average change in RMSE between these two σw values was found to be 5.5%.

Discussion

The residual plots over time for the fitted ODE model shown in Fig. 3 suggest that there is correlation among the residuals for most of the occasions and that the errors are not arising due to randomness. These correlated residuals should be accounted for by modelling with SDEs.

Thus, we used an SDE pharmacokinetic model of LCIG for patients with Parkinson’s disease in our study. The parameter estimated in this study is bioavailability, which is the proportion of levodopa that enters the circulation when introduced into the body. When all 20 occasions are used, both models estimated similar values for bioavailability at around 75%, but there are too few patients for reliable population parameter estimation. In the study by Westin et al. [6], there was an additional 36 occasions with generally stable LCIG concentrations – that would increase the computation time significantly for the cross-validation performed, which is a key objective of our study to compare model fits. The variability of bioavailability is estimated to be too small and close to 0, which is in accordance with Westin et al. [6], i.e., that a fixed bioavailability parameter is sufficient to fit individual models. This parameter was also fixed to 0 in a study by Thomas et al. [21], where they studied individual dose–response models for levodopa infusion dose optimisation.

When a stochastic term is introduced in Eq. 9, the diffusion scaling term, σw, is estimated to be 0.099. This is significantly different from zero, implying that the intra-individual variability has been separated into measurement and system noises.

A cross-validation study was also performed for both ODE and SDE models to check how well models predict data, which had not been performed previously. From Fig. 4, it is seen that introducing stochasticity does decrease the average RMSE. When σw is fixed to 0 (which is the ODE model), the average RMSE (0.313) is found to be the highest; when σw is fixed to 0.9 (which is an SDE model), the average RMSE (0.297) is found to be lowest. The binomial test suggests that these RMSEs are significantly different, so the SDE model does fit the data better.

Both the ODE and SDE models are reasonably effective at fitting observed data. It was found that the mean percentage change in RMSE was 5.5% between the two aforementioned σw values. Even though visual checks do not provide a distinguishable difference between the ODE and SDE models, the cross-validation study suggests that the SDE model is 5.5% better at predicting the LCIG plasma levodopa concentration in the system, based on dose information and individualised parameters.

The maximum likelihood estimate for σw was estimated to be 0.099, however, where the average RMSE was not the lowest according to the cross-validation study. This could be because the true value of the global minimum has not been found, as suggested by Miller et al. [14]. A cross-validation study may therefore be useful in investigating an appropriate value for σw in an SDE model. Other options include using different optimisation algorithms, or using a Bayesian approach for pharmacokinetic/pharmacodynamic modelling, which were not explored in our study.

The cross-validation method that we used in the present study could be applied in other mixed-effects models, not limited to clinical studies, where there are time series data from multiple occasions for each subject. By leaving out all but one occasion for a subject in the training set, random-effect estimates can be extracted for that particular subject, which are then fixed in the test set (the occasions left out).

A limitation of our study is that almost all the parameters were fixed to population values found in the literature (while allowing for individual variations which have been estimated). This was done to reduce the computation time for the overall analysis, as cross-validation was performed for 11 different values of σw, and there were 34 parameter estimation computations at each of those values.

A model was tested where the σw term was included in all three equations (8–10) and they were estimated. The estimations yielded equal values for σw, but the computation demands for this model were rather intensive. Further analyses were performed to investigate the necessity of adding a σw term in Eq. 8 by simulating stochastic noise in the first compartment. This gave very similar predictions of the concentration in the second compartment, regardless of whether the stochastic noise was assumed to be in Eq. 8 or 9. An additional parameter, TABS, was also estimated along with bioavailability to compare their estimates with the ODE model and to investigate the value of the σw term. The parameter estimates were similar and the results were identical to the simpler model shown in this study.

It is to be noted, however, that infusion data are generally stable and therefore may not show much volatility. A different dataset, such as one with oral dosage information, could yield different results. It would also be interesting to study whether stochasticity may be included in the pharmacodynamic effect equation, where the observations are more volatile, i.e., the Wiener component could be introduced into the effect equation of the pharmacokinetic/pharmacodynamic model of LCIG levodopa as described by Westin et al. [6]. It could also be investigated whether an effect prediction may be fitted directly from dose information in an SDE setting.

In conclusion, we investigated the pharmacokinetic model of levodopa infusion using SDEs. Correlation between the residuals was observed, so a stochastic Wiener component was added to the system to account for auto-correlated errors. The SDE model provided better fits and was 5.5% better, in terms of lower RMSE, at predicting LCIG levodopa pharmacokinetics based on dose information.

Stochastic differential equations modelling can provide better predictions for plasma levodopa concentration in patients with Parkinson’s disease.

Cross-validation can be used to investigate the value of the diffusion scale parameter that provides the best fit in a stochastic differential equations model.