Separate Vs Simultaneous Fitting

When we want to model both the time course of the drug in the plasma and the effect it has on a biomarker, there are two options open to us.

We can write one simultaneous model where the concentration of a drug in one compartment directly affects the biomarker in another.

Alternatively, we can write one model for the PK part and then a separate model for the PD part, using estimates from the first part to fix parameters in the PD model.

Both strategies have advantages and disadvantages in different situations.

In most cases, a separate fit is used because it takes less computational time and the results are similar (http://www.ncbi.nlm.nih.gov/pubmed/15000421).

However, when the PK data is sparse (few measurements per individual) and the PD data is rich (more measurements per individual), the extra information provided by the PD data may lead to a more accurate model if a simultaneous fit is used. In a real modelling situation we would not know the generating values, in contrast to the following examples, which means we do not know which method is more accurate. If there are are more PD observations than PK observations it is worth considering a simultaneous fit.

Here are some examples of PKPD model where models are fitted simultaneously, then both PK and PD models are fit separately.

For each model, we then take away some of the PK data to see how this affects the simultaneous and separate fits.

Simple PKPD Model

Simultaneous fit with same number of observations

The PK part is a simple one compartment model.

The concentration in the central compartment affects the rate at which the biomarker enters the body. A larger concentration will lead to an increased rate of entry into the body.

There are 40 individuals and each individual has 10 observations of PK and PD measurements.

The DERIVATIVES block defines the equations for the PK and PD compartments:

DERIVATIVES: |

    d[CENTRAL] = @bolus{amt: c[DOSE]} - s[CENTRAL]*m[CL]/m[V]
    conc = s[CENTRAL]/m[V]
    d[BIO_CENTRAL] = m[MARKER] * conc - m[BIO_K] * s[BIO_CENTRAL]

The PREDICTIONS block has two parts to it because the are two measurements of interest; the PK measurements (defined as DV_CENTRAL in this script) and the PD measurements (defined as BIO_CENTRAL in this script):

PREDICTIONS: |

    predlabel[DV_CENTRAL] = 'Drug conc'
    pred[DV_CENTRAL] = s[CENTRAL]/m[V]
    dv_central_var = (pred[DV_CENTRAL]*m[PNOISEPK])**2 + m[ANOISEPK]**2
    c[DV_CENTRAL] ~ norm(pred[DV_CENTRAL], dv_central_var)

    predlabel[BIO_CENTRAL] = 'Biomarker'
    pred[BIO_CENTRAL] = s[BIO_CENTRAL]
    bio_central_var = (pred[BIO_CENTRAL]*m[PNOISEPD])**2 + m[ANOISEPD]**2
    c[BIO_CENTRAL] ~ norm(pred[BIO_CENTRAL], bio_central_var)

The compartment diagram is:

PK fit with same number of observations

We used the following script in the R programme (https://www.r-project.org/) to remove the PD data from the synthetic_data.csv file:

orig<- read.csv("synthetic_data.csv")
orig<- orig[orig$CMT_LABEL == "CENTRAL", ]
write.csv(orig, file = "PK_data.csv", quote = F, row.names=F)

The DERIVATIVES block now only has one compartment:

DERIVATIVES: |

    d[CENTRAL] = @bolus{amt: c[DOSE]} - s[CENTRAL]*m[CL]/m[V]

And the PREDICTIONS block has one part to it:

PREDICTIONS: |

    predlabel[DV_CENTRAL] = 'Drug conc'
    pred[DV_CENTRAL] = s[CENTRAL]/m[V]
    dv_central_var = (pred[DV_CENTRAL]*m[PNOISEPK])**2 + m[ANOISEPK]**2
    c[DV_CENTRAL] ~ norm(pred[DV_CENTRAL], dv_central_var)

PD fit with same number of observations

We used the following R script to remove the PK data from the synthetic_data.csv file and add the estimated PK parameters from the px_params.csv file:

orig<- read.csv("synthetic_data.csv")
pk<- read.csv("px_params.csv")
orig<- orig[orig$CMT_LABEL == "BIO_CENTRAL"|orig$TYPE == "reset"|orig$TYPE == "dose", ]
orig$CL<- pk$CL
orig$V<- pk$V
write.csv(orig, file = "PD_data.csv", quote = F, row.names=F)

This is reflected in the script file by using the values from the data file (c[CL] and c[V]) instead of predicting the parameters m[CL] and m[V]:

DERIVATIVES: |

    d[CENTRAL] = @bolus{amt: c[DOSE]} - s[CENTRAL]*c[CL]/c[V]
    conc = s[CENTRAL]/c[V]
    d[BIO_CENTRAL] = m[MARKER] * conc - m[BIO_K] * s[BIO_CENTRAL]

The PREDICTIONS block has only one part because only the PD measurements are being fitted:

PREDICTIONS: |

    predlabel[BIO_CENTRAL] = 'Biomarker'
    pred[BIO_CENTRAL] = s[BIO_CENTRAL]
    bio_central_var = (pred[BIO_CENTRAL]*m[PNOISEPD])**2 + m[ANOISEPD]**2
    c[BIO_CENTRAL] ~ norm(pred[BIO_CENTRAL], bio_central_var)

The differences in each estimated parameter for each different method are summarised in the table below:

Parameter	Generation Value	Simultaneous estimate	Percent diff from original	Separate estimate	Percent diff from original
CL	2	2.01	0.5	2.00	0
V	25	25.58	2.3	25.49	2.0
MARKER	0.3	0.284	-5.3	0.284	-5.3
BIO_K	5	5.16	3.2	5.14	2.8
PNOISEPK	0.06	0.033	-45	0.049	-18.3
ANOISEPK	0.5	0.473	-5.4	0.406	-18.8
PNOISEPD	0.03	0.030	0	0.033	10
ANOISEPD	0.5	0.560	12	0.568	13.6
CL_isv	0.02	0.020	0	0.020	0
V_isv	0.04	0.037	-7.5	0.039	-2.5
MARKER_isv	0.02	0.020	0	0.026	30
BIO_K_isv	0.05	0.050	0	0.056	12.4

There is very little difference is the parameters estimated by the two methods.

Simultaneous fit with sparse PK data and rich PD data

orig<- read.csv("synthetic_data.csv")
orig<- orig[orig$CMT_LABEL == "BIO_CENTRAL"|orig$TYPE == "reset"|orig$TYPE == "dose"| orig$TIME == 4.9| orig$TIME == 24.4, ] # removes all but 2 PK data
write.csv(orig, file = "PKPD_data.csv", quote = F, row.names=F)

PK fit with sparse PK data

We used the following R script to remove the PD data from the PKPD_data.csv file:

orig<- read.csv("PKPD_data.csv")
orig<- orig[orig$CMT_LABEL == "CENTRAL", ]
write.csv(orig, file = "PK_data.csv", quote = F, row.names=F)

PD fit with sparse PK data and rich PD data

We used the following R script to remove the PK data from the PKPD_data.csv file and add the estimated PK parameters from the px_params.csv file:

orig<- read.csv("PKPD_data.csv")
pk<- read.csv("px_params.csv")
orig<- orig[orig$CMT_LABEL == "BIO_CENTRAL"|orig$TYPE == "reset"|orig$TYPE == "dose", ]
ind<- match(orig$ID, pk$ID)
orig$CL<- pk$CL[ind]
orig$V<- pk$V[ind]

The differences in each estimated parameter for the simultaneous and separate methods with the sparse PK and rich PD data are summarised in the table below:

Parameter	Generation Value	Simultaneous estimate	Percent diff from original	Separate estimate	Precent diff from original
CL	2	2.01	0.5	2.02	1.0
V	25	25.61	2.4	25.56	2.2
MARKER	0.3	0.284	-5.3	0.287	-4.3
BIO_K	5	5.16	3.2	5.20	4
PNOISEPK	0.06	0.045	-25	0.058	-3.3
ANOISEPK	0.5	0.522	4.4	0.892	78
PNOISEPD	0.03	0.030	0	0.045	50
ANOISEPD	0.5	0.560	12	0.615	23
CL_isv	0.02	0.025	25	0.018	-10
V_isv	0.04	0.036	-10	0.032	-20
MARKER_isv	0.02	0.018	5	0.059	195
BIO_K_isv	0.05	0.049	4	0.072	44

There is very little difference in the PK and PD parameter estimates between the two methods, but there are large differences in the estimates of noise and inter-individual variation. The sparse PK data has led to larger estimates of the noise in the model when separate fitting is used in contrast to the simultaneous fit. The PD parameters also have much larger estimates of inter-individual variation.

In the situation of having sparse PK data and rich PD data, there is an advantage to a simultaneous fit compared to a separate fit using this model.

Emax PKPD Model

Emax Simultaneous fit with same number of observations

In this case, the PK part is a three compartment model.

The concentration in the central compartment affects the rate at which a biomarker leaves the body using an emax equation.

There are 40 individuals and each individual has 10 observations of PK and PD measurements.

Many of the parameters used in the DERIVATIVES block, such as the flows between compartments (eg. m[K12]) or the flow into the PD compartment (m[KIN]), are defined in the MODEL_PARAMS section to keep the equations in the DERIVATIVES block shorter:

MODEL_PARAMS: |

    m[CL] = f[CL] * exp(r[CL])
    m[V1] = f[V1] * exp(r[V1])
    m[Q2] = f[Q2] * exp (r[Q2])
    m[V2] = f[V2] * exp(r[V2])
    m[Q3] = f[Q3] * exp (r[Q3])
    m[V3] = f[V3] * exp(r[V3])
    m[K] = m[CL]/m[V1]
    m[K12] = m[Q2]/m[V1]
    m[K13] =  m[Q3]/m[V1]
    m[K21] =  m[Q2]/m[V2]
    m[K31] =  m[Q3]/m[V3]
    m[BASE] = f[BASE] * exp(r[BASE])
    m[KOUT] = f[KOUT] * exp(r[KOUT])
    m[KIN] = m[BASE]*m[KOUT]
    m[C50] = f[C50] * exp(r[C50])
    m[EMAX] = f[EMAX] * exp(r[EMAX])
    m[PNOISEPK] = f[PNOISEPK]
    m[ANOISEPK] = f[ANOISEPK]
    m[PNOISEPD] = f[PNOISEPD]
    m[ANOISEPD] = f[ANOISEPD]

The DERIVATIVES block defines the equations for the PK and PD compartments:

DERIVATIVES: |

    d[CENTRAL] = @bolus{amt: c[DOSE]} - m[K12]*s[CENTRAL] + m[K21]*s[PERI1] - m[K13]*s[CENTRAL] + m[K31]* s[PERI2] - m[K]*s[CENTRAL]
    conc = s[CENTRAL]/m[V1]
    d[PERI1] = m[K12]*s[CENTRAL] - m[K21]*s[PERI1]
    d[PERI2] = m[K13]*s[CENTRAL] - m[K31]*s[PERI2]
    eff = m[EMAX]* conc/(m[C50]+ conc)
    d[BIO_CENTRAL] = m[KIN] - m[KOUT]*s[BIO_CENTRAL]* eff

The PREDICTIONS is the same as in Simple PKPD Model with two parts to it because there are two measurements of interest; the PK and the PD.

The compartment diagram is:

Emax PK fit with same number of observations

We used the same R script as in PK fit with same number of observations to remove the PD data from the synthetic_data.csv file.

The DERIVATIVES block now has three compartments:

DERIVATIVES: |

    d[CENTRAL] = @bolus{amt: c[DOSE]} - m[K12]*s[CENTRAL] + m[K21]*s[PERI1] - m[K13]*s[CENTRAL] + m[K31]* s[PERI2] - m[K]*s[CENTRAL]
    d[PERI1] = m[K12]*s[CENTRAL] - m[K21]*s[PERI1]
    d[PERI2] = m[K13]*s[CENTRAL] - m[K31]*s[PERI2]

Emax PD fit with same number of observations

We used the following R script to remove the PK data from the synthetic_data.csv file and add the estimated PK parameters from the px_params.csv file:

orig<- read.csv("synthetic_data.csv")
pk<- read.csv("px_params.csv")

orig<- orig[orig$CMT_LABEL == "BIO_CENTRAL"|orig$TYPE == "reset"|orig$TYPE == "dose", ]
orig$CL<- pk$CL
orig$V1<- pk$V1
orig$Q2<- pk$Q2
orig$V2<- pk$V2
orig$Q3<- pk$Q3
orig$V3<- pk$V3

write.csv(orig, file = "PD_data.csv", quote = F, row.names=F)

Some of the parameters used in the PK section of the DERIVATIVES block are defined in the MODEL_PARAMS section using fixed values for each individual from the data file instead of estimating them as was the case in the simultaneous fit:

MODEL_PARAMS: |

    m[K] = c[CL]/c[V1]
    m[K12] = c[Q2]/c[V1]
    m[K13] =  c[Q3]/c[V1]
    m[K21] =  c[Q2]/c[V2]
    m[K31] =  c[Q3]/c[V3]
    m[BASE] = f[BASE] * exp(r[BASE])
    m[KOUT] = f[KOUT] * exp(r[KOUT])
    m[KIN] = m[BASE]*m[KOUT]
    m[C50] = f[C50] * exp(r[C50])
    m[EMAX] = f[EMAX] * exp(r[EMAX])
    m[PNOISEPD] = f[PNOISEPD]
    m[ANOISEPD] = f[ANOISEPD]

The DERIVATIVES section is very similar to the simultaneous fit, except the concentration is defined by dividing the amount in the central compartment by the fixed value of volume of distribution for each individual:

DERIVATIVES: |

    d[CENTRAL] = @bolus{amt: c[DOSE]} - m[K12]*s[CENTRAL] + m[K21]*s[PERI1] - m[K13]*s[CENTRAL] + m[K31]* s[PERI2] - m[K]*s[CENTRAL]
    conc = s[CENTRAL]/c[V1]
    d[PERI1] = m[K12]*s[CENTRAL] - m[K21]*s[PERI1]
    d[PERI2] = m[K13]*s[CENTRAL] - m[K31]*s[PERI2]
    eff = m[EMAX]* conc/(m[C50]+ conc)
    d[BIO_CENTRAL] = m[KIN] - m[KOUT]*s[BIO_CENTRAL]* eff

The PREDICTIONS block has only one part because only the PD measurements are being fitted:

PREDICTIONS: |

    predlabel[BIO_CENTRAL] = 'Biomarker'
    pred[BIO_CENTRAL] = s[BIO_CENTRAL]
    bio_central_var = (pred[BIO_CENTRAL]*m[PNOISEPD])**2 + m[ANOISEPD]**2
    c[BIO_CENTRAL] ~ norm(pred[BIO_CENTRAL], bio_central_var)

The differences in each estimated parameter for each different method are summarised in the table below:

Parameter	Generation Value	Simultaneous estimate	Percent diff from original	Separate estimate	Percent diff from original
CL	1.5	1.69	12.7	1.69	12.7
V1	15	15.26	1.7	15.25	1.7
Q2	3	2.74	-8.7	2.75	-8.3
V2	10	11.52	15.2	11.77	17.7
Q3	0.5	0.23	-54	0.21	-58
V3	25	14.87	-40.5	15.02	-39.9
BASE	100	111.53	11.5	107.96	8.0
KOUT	0.2	0.18	-10	0.18	-10
C50	20	17.5	-5.3	0.284	-12.5
EMAX	50	50.3	0.6	51.8	3.6
PNOISEPK	0.05	0.044	-12	0.046	-8
ANOISEPK	0.5	0.471	-5.8	0.463	-7.4
PNOISEPD	0.05	0.049	-2	0.052	4
ANOISEPD	1	1.01	1	1.01	1
CL_isv	0.05	0.038	-24	0.039	-22
V1_isv	0.02	0.025	25	0.025	25
Q2_isv	0.05	0.017	-66 `	0.0089	-82.2
V2_isv	0.02	0.017	-15	0.016	-20
Q3_isv	0.01	0.034	240	0.047	370
V3_isv	0.01	0.035	250	0.031	210
BASE_isv	0.05	0.026	-48	0.027	-46
KOUT_isv	0.02	0.030	50	0.035	75
C50_isv	0.03	0.099	230	0.12	300
EMAX_isv	0.01	0.014	40	0.016	60

There is very little difference is the parameters estimated by the two methods.

Emax Simultaneous fit with sparse PK data and rich PD data

However, the input data file used for the model has been altered so that there are only 2 PK measurements per individual whilst retaining the 10 PD measurements per individual.

Emax PK fit with sparse PK data

We used the R script in PK fit with sparse PK data to remove the PD data from the PKPD_data.csv file.

PD fit with sparse PK data and rich PD data

We used the following R script to remove the PK data from the PKPD_data.csv file:-

orig<- read.csv("PKPD_data.csv")
pk<- read.csv("px_params.csv")

orig<- orig[orig$CMT_LABEL == "BIO_CENTRAL"|orig$TYPE == "reset"|orig$TYPE == "dose", ]
ind<- match(orig$ID, pk$ID)
orig$CL<- pk$CL[ind]
orig$V1<- pk$V1[ind]
orig$Q2<- pk$Q2[ind]
orig$V2<- pk$V2[ind]
orig$Q3<- pk$Q3[ind]
orig$V3<- pk$V3[ind]

write.csv(orig, file = "PD_data.csv", quote = F, row.names=F)

The differences in each estimated parameter for the simultaneous and separate methods with the sparse PK and rich PD data are summarised in the table below:

Parameter	Generation Value	Simultaneous estimate	Percent diff from original	Separate estimate	Percent diff from original
CL	1.5	1.76	17.4	1.62	14.4
V1	15	20.1	34	19	26.4
Q2	3	2.09	-30.5	2.75	-34.6
V2	10	8.10	-19	10.7	6.8
Q3	0.5	0.20	-59.5	0.21	-58
V3	25	15.5	-38	15.3	-38.7
BASE	100	114.63	14.6	103.7	3.7
KOUT	0.2	0.17	-13.8	0.19	-4.7
C50	20	15.4	-38	15.3	-38.7
EMAX	50	47.5	-5	53.8	7.6
PNOISEPK	0.05	0.024	-52.8	0.019	-62.8
ANOISEPK	0.5	0.266	-46.8	0.415	-17
PNOISEPD	0.05	0.056	11.4	0.053	5.6
ANOISEPD	1	0.99	-0.7	1	0.4
CL_isv	0.05	0.043	-13.4	0.052	4.2
V1_isv	0.02	0.018	-8	0.035	75
Q2_isv	0.05	0.077	54.2	0.034	-31.2
V2_isv	0.02	0.025	24	0.031	55.5
Q3_isv	0.01	0.025	149	0.031	208
V3_isv	0.01	0.028	176	0.03	201
BASE_isv	0.05	0.031	-37.4	0.028	-43.6
KOUT_isv	0.02	0.028	37.5	0.037	83
C50_isv	0.03	0.059	95	0.128	327
EMAX_isv	0.01	0.012	20	0.018	77

In contrast to Simple PKPD Model most parameters have estimates closer to the generating value when they are estimated separately as opposed to simultaneously. The inter-individual variance is generally much larger when the separate fitting is used, especially for the PD parameters. In this case, there does not appear to be an advantage to using simultaneous fitting over separate fitting despite there being more information in the PD measurements.

Absorption PKPD Model

The PK part is a two compartment model with first order absorption.

The concentration in the central compartment affects the rate at which the biomarker enters the body. A larger concentration will lead to an increased rate of entry into the body.

There are 30 individuals and each individual has 20 observations of PK and PD measurements.

The DERIVATIVES block defines the equations for the PK and PD compartments:

DERIVATIVES: |
    d[DEPOT] = @bolus{amt: c[DOSE]} - m[KA]*s[DEPOT]
    d[CENTRAL] = m[KA]*s[DEPOT] - s[CENTRAL]*m[Q]/m[V1] + s[PERI]*m[Q]/m[V2] - s[CENTRAL]*m[CL]/m[V1]
    conc = s[CENTRAL]/m[V1]
    d[PERI] = s[CENTRAL]*m[Q]/m[V1] - s[PERI]*m[Q]/m[V2]
    d[BIO_CENTRAL] = m[MARKER] * conc - m[BIO_K] * s[BIO_CENTRAL]

The PREDICTIONS block has two parts to it because the are two measurements of interest; the PK measurements (defined as DV_CENTRAL in this script) and the PD measurements (defined as BIO_CENTRAL in this script):

PREDICTIONS: |
    predlabel[DV_CENTRAL] = 'Drug conc'
    pred[DV_CENTRAL] = s[CENTRAL]/m[V1]
    dv_central_var = (pred[DV_CENTRAL]*m[PNOISEPK])**2 + m[ANOISEPK]**2
    c[DV_CENTRAL] ~ norm(pred[DV_CENTRAL], dv_central_var)

    predlabel[BIO_CENTRAL] = 'Biomarker'
    pred[BIO_CENTRAL] = s[BIO_CENTRAL]
    bio_central_var = (pred[BIO_CENTRAL]*m[PNOISEPD])**2 + m[ANOISEPD]**2
    c[BIO_CENTRAL] ~ norm(pred[BIO_CENTRAL], bio_central_var)

The compartment diagram is:

Absorption PK fit with same number of observations

orig<- read.csv("synthetic_data.csv")
orig<- orig[orig$CMT_LABEL == "CENTRAL"|orig$CMT_LABEL == "DEPOT", ]
write.csv(orig, file = "PK_data.csv", quote = F, row.names=F)

The DERIVATIVES block now only has two compartments and an absorption compartment:

DERIVATIVES: |

    d[DEPOT] = @bolus{amt: c[DOSE]} - m[KA]*s[DEPOT]
    d[CENTRAL] = m[KA]*s[DEPOT] - s[CENTRAL]*m[Q]/m[V1] + s[PERI]*m[Q]/m[V2] - s[CENTRAL]*m[CL]/m[V1]
    d[PERI] = s[CENTRAL]*m[Q]/m[V1] - s[PERI]*m[Q]/m[V2]

And the PREDICTIONS block has one part to it:

PREDICTIONS: |

    predlabel[DV_CENTRAL] = 'Drug conc'
    pred[DV_CENTRAL] = s[CENTRAL]/m[V1]
    dv_central_var = (pred[DV_CENTRAL]*m[PNOISEPK])**2 + m[ANOISEPK]**2
    c[DV_CENTRAL] ~ norm(pred[DV_CENTRAL], dv_central_var)

Absorption PD fit with same number of observations

orig<- read.csv("synthetic_data.csv")
pk<- read.csv("px_params.csv")
orig<- orig[orig$CMT_LABEL == "BIO_CENTRAL"|orig$TYPE == "reset"|orig$TYPE == "dose", ]
orig$KA<- pk$KA
orig$CL<- pk$CL
orig$V1<- pk$V1
orig$Q<- pk$Q
orig$V2<- pk$V2
write.csv(orig, file = "PD_data.csv", quote = F, row.names=F)

This is achieved in the script file by using the values from the data file (c[KA], c[CL], c[V1], c[Q] and c[V2]) instead of predicting these parameters:

DERIVATIVES: |

d[DEPOT] = @bolus{amt: c[DOSE]} - c[KA]*s[DEPOT]
d[CENTRAL] = c[KA]*s[DEPOT] - s[CENTRAL]*c[Q]/c[V1] + s[PERI]*c[Q]/c[V2] - s[CENTRAL]*c[CL]/c[V1]
conc = s[CENTRAL]/c[V1]
d[PERI] = s[CENTRAL]*c[Q]/c[V1] - s[PERI]*c[Q]/c[V2]
d[BIO_CENTRAL] = m[MARKER] * conc - m[BIO_K] * s[BIO_CENTRAL]

The PREDICTIONS block has only one part because only the PD measurements are being fitted:

PREDICTIONS: |

    predlabel[BIO_CENTRAL] = 'Biomarker'
    pred[BIO_CENTRAL] = s[BIO_CENTRAL]
    bio_central_var = (pred[BIO_CENTRAL]*m[PNOISEPD])**2 + m[ANOISEPD]**2
    c[BIO_CENTRAL] ~ norm(pred[BIO_CENTRAL], bio_central_var)

The differences in each estimated parameter for each different method are summarised in the table below:

Parameter	Generation Value	Simultaneous estimate	Percent diff from original	Separate estimate	Percent diff from original
KA	1.5	1.50	0.03	1.40	-6.9
CL	2	1.92	-4.1	1.88	-5.9
V1	25	26.34	5.3	26.16	4.7
Q	1	1.00	0.3	1.04	4.4
V2	15	16.36	9.1	18.03	20.2
BIO_K	0.3	0.31	3.4	0.31	4.3
MARKER	5	5.32	6.4	5.38	7.7
PNOISEPK	0.02	0.021	4	0.024	21.5
ANOISEPK	0.5	0.508	1.5	0.502	0.34
PNOISEPD	0.03	0.03	0	0.036	19.3
ANOISEPD	0.5	0.558	11.6	0.649	29.7
KA_isv	0.03	0.033	7.7	0.017	-43.7
CL_isv	0.02	0.019	-7.5	0.025	22.5
V1_isv	0.04	0.042	4.3	0.047	18.3
Q_isv	0.01	0.029	192	0.012	15
V2_isv	0.05	0.038	-24.4	0.132	163
BIO_K_isv	0.02	0.032	57.5	0.034	69
MARKER_isv	0.05	0.059	17.8	0.060	19.4

In this case, the estimates from the simultaneous method are generally closer to the generating values than when the PK and PD parameters are estimated separately.

Absorption Simultaneous fit with sparse PK data and rich PD data

orig<- read.csv("synthetic_data.csv")
orig<- orig[orig$CMT_LABEL == "BIO_CENTRAL"|orig$TYPE == "reset"|orig$TYPE == "dose"| orig$TIME == 2.95| orig$TIME == 12.70| orig$TIME == 32.20, ] # removes all but 3 PK data
write.csv(orig, file = "PKPD_data.csv", quote = F, row.names=F)

Absorption PK fit with sparse PK data

orig<- read.csv("PKPD_data.csv")
pk<- read.csv("px_params.csv")
orig<- orig[orig$CMT_LABEL == "BIO_CENTRAL"|orig$TYPE == "reset"|orig$TYPE == "dose", ]
ind<- match(orig$ID,pk$ID)
orig$KA<- pk$KA[ind]
orig$CL<- pk$CL[ind]
orig$V1<- pk$V1[ind]
orig$Q<- pk$Q[ind]
orig$V2<- pk$V2[ind]

The differences in each estimated parameter for the simultaneous and separate methods with the sparse PK and rich PD data are summarised in the table below:

Parameter	Generation Value	Simultaneous estimate	Percent diff from original	Separate estimate	Percent diff from original
KA	1.5	1.54	2.9	1.62	8.3
CL	2	1.92	-4.1	1.80	-9.8
V1	25	26.34	5.3	26.96	7.8
Q	1	1.00	0.3	1.01	0.8
V2	15	16.36	9.1	22.62	50.8
BIO_K	0.3	0.31	3.4	0.31	3.3
MARKER	5	5.32	6.4	5.31	6.2
PNOISEPK	0.02	0.021	4	0.045	125.5
ANOISEPK	0.5	0.508	1.5	0.001	-99.8
PNOISEPD	0.03	0.03	0	0.040	33.7
ANOISEPD	0.5	0.558	11.6	2.38	376.9
KA_isv	0.03	0.033	7.7	0.054	78.3
CL_isv	0.02	0.019	-7.5	0.030	48
V1_isv	0.04	0.042	4.3	0.041	2.3
Q_isv	0.01	0.029	192	0.02	100
V2_isv	0.05	0.038	-24.4	0.159	218
BIO_K_isv	0.02	0.032	57.5	0.045	124
MARKER_isv	0.05	0.059	17.8	0.057	13.8

Unsurprisingly, the simultaneous estimates are again closer to the generating values than the separate estimates. One of the parameters (ANOISEPK) has been pushed downwards to its minimum value of 0.001 in the separate fit, which gives us an indication that the model is not a good fit to the data. In this case, it may be because there are not enough PK observations for an acceptable fit.

Absorption PD fit with sparse PK data and rich PD data

orig<- read.csv("PKPD_data.csv")
pk<- read.csv("px_params.csv")
orig<- orig[orig$CMT_LABEL == "BIO_CENTRAL"|orig$TYPE == "reset"|orig$TYPE == "dose", ]
ind<- match(orig$ID,pk$ID)
orig$KA<- pk$KA[ind]
orig$CL<- pk$CL[ind]
orig$V1<- pk$V1[ind]
orig$Q<- pk$Q[ind]
orig$V2<- pk$V2[ind]

The differences in each estimated parameter for the simultaneous and separate methods with the sparse PK and rich PD data are summarised in the table below:

Parameter	Generation Value	Simultaneous estimate	Percent diff from original	Separate estimate	Percent diff from original
KA	1.5	1.54	2.9	1.62	8.3
CL	2	1.92	-4.1	1.80	-9.8
V1	25	26.34	5.3	26.96	7.8
Q	1	1.00	0.3	1.01	0.8
V2	15	16.36	9.1	22.62	50.8
BIO_K	0.3	0.31	3.4	0.31	3.3
MARKER	5	5.32	6.4	5.31	6.2
PNOISEPK	0.02	0.021	4	0.045	125.5
ANOISEPK	0.5	0.508	1.5	0.001	-99.8
PNOISEPD	0.03	0.03	0	0.040	33.7
ANOISEPD	0.5	0.558	11.6	2.38	376.9
KA_isv	0.03	0.033	7.7	0.054	78.3
CL_isv	0.02	0.019	-7.5	0.030	48
V1_isv	0.04	0.042	4.3	0.041	2.3
Q_isv	0.01	0.029	192	0.02	100
V2_isv	0.05	0.038	-24.4	0.159	218
BIO_K_isv	0.02	0.032	57.5	0.045	124
MARKER_isv	0.05	0.059	17.8	0.057	13.8

Unsurprisingly, the simultaneous estimates are again closer to the generating values than the separate estimates. One of the parameters (ANOISEPK) has been pushed downwards to its minimum value of 0.001 in the separate fit, which gives us an indication that the model is not a good fit to the data. In this case, it may be because there are not enough PK observations for an acceptable fit.