using DeepPumas
using CairoMakie
using StableRNGs
using PumasPlots
set_mlp_backend(:staticflux)
Mixed Effects Pharmacodynamics Identification using DeepNLME
1 Learning objectives
By the end of this tutorial, you will be able to use neural-embedded nonlinear mixed effects (DeepNLME) modeling to learn and improve existing dynamics models in the presence of random effects.
2 Setup
First, we need to load the required packages. We also set the backend for the MLPDomain to use :staticflux.
3 Synthetic data generation
We will generate synthetic data from the indirect response model (IDR) model defined below.
datamodel = @model begin
@param begin
tvKa ∈ RealDomain()
tvCL ∈ RealDomain()
tvVc ∈ RealDomain()
tvSmax ∈ RealDomain()
tvn ∈ RealDomain()
tvSC50 ∈ RealDomain()
tvKout ∈ RealDomain()
tvKin ∈ RealDomain()
Ω ∈ PDiagDomain(5)
σ ∈ RealDomain()
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
Smax = tvSmax * exp(η[1])
SC50 = tvSC50 * exp(η[2])
Ka = tvKa * exp(η[3])
Vc = tvVc * exp(η[4])
Kout = tvKout * exp(η[5])
Kin = tvKin
CL = tvCL
n = tvn
end
@init begin
Response = Kin / Kout
end
@vars begin
cp := max(Central / Vc, 0.0)
EFF := Smax * cp^n / (SC50^n + cp^n)
end
@dynamics begin
Depot' = -Ka * Depot
Central' = Ka * Depot - (CL / Vc) * Central
Response' = Kin * (1 + EFF) - Kout * Response
end
@derived begin
Outcome ~ @. Normal(Response, σ)
end
endThe following parameters are used to generate the synthetic data.
p_data = (;
tvKa = 0.5,
tvCL = 1.0,
tvVc = 1.0,
tvSmax = 2.9,
tvn = 1.5,
tvSC50 = 0.05,
tvKout = 2.2,
tvKin = 0.8,
Ω = Diagonal(fill(0.1, 5)),
σ = 0.1,
)To generate the synthetic data, we will use the simobs function. The obstimes argument specifies the time points at which observations are made. We assume all subjects are observed at times 0, 1, 2, …, 24 hours. Each subject receives two doses of 1 unit. The first dose is given at baseline and the second dose is at a random time centered around hour 5. A total of 22 subjects are generated, with 20 subjects used for training and 12 for testing.
obstimes = 0:24
ntrain = 20
ntest = 12
pop = map(1:ntrain+ntest) do i
rng = StableRNG(i)
dose_1 = DosageRegimen(1.0)
dose_2 = DosageRegimen(1.0; time = rand(rng, Gamma(40, 5 / 40)))
sim = simobs(
datamodel,
Subject(; id = i, events = DosageRegimen(dose_1, dose_2)),
p_data;
obstimes,
rng,
)
Subject(sim)
end
trainpop = pop[1:ntrain]
testpop = pop[(ntrain+1):end]We then visualize the synthetic data and the predictions of the data-generating model. The specified obstimes with a smaller time step is just to get a smooth plot.
pred_datamodel = predict(datamodel, testpop, p_data; obstimes = 0:0.1:24);
plotgrid(pred_datamodel)
4 Neural-embedded NLME modeling
To model the synthetic data, we will use the neural-embedded nonlinear mixed effects (DeepNLME) model shown below. This model includes a multi-layer perceptron (MLP) to capture the nonlinear relationship between the drug concentration and the response. The MLP is defined in the @param block using the MLPDomain function. The MLP has 5 inputs (2 state variables + 3 random effects) and a single output.
model = @model begin
@param begin
NN ∈ MLPDomain(5, 7, 7, (1, identity); reg = L2(1.0))
tvKa ∈ RealDomain(; lower = 0)
tvCL ∈ RealDomain(; lower = 0)
tvVc ∈ RealDomain(; lower = 0)
tvR₀ ∈ RealDomain(; lower = 0)
ωR₀ ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ ∈ RealDomain(; lower = 0)
end
@random begin
η ~ MvNormal(Ω)
η_nn ~ MvNormal(0.01 * I(3))
end
@pre begin
Ka = tvKa * exp(η[1])
Vc = tvVc * exp(η[2])
CL = tvCL
Response₀ = tvR₀ * exp(10 * ωR₀ * η_nn[1])
iNN = fix(NN, η_nn)
end
@init begin
Response = Response₀
end
@dynamics begin
Depot' = -Ka * Depot
Central' = Ka * Depot - (CL / Vc) * Central
Response' = iNN(Central / Vc, Response)[1]
end
@derived begin
Outcome ~ @. Normal(Response, σ)
end
endThe model above has 5 random effects, 2 in η and 3 in η_nn. The two random effects in η are used to model the between-subject variability of the pharmacokinetic (PK) parameters. The remaining 3 random effects in η_nn are used as input to the MLP. However, η_nn[1] is also used to define the initial value of the response variable Response using Response₀ = tvR₀ * exp(10 * ωR₀ * η_nn[1]).
Note that η_nn[1] is a normally distributed random variable with mean 0 and standard deviation 1. The term 10 * ωR₀ * η_nn[1] is therefore another normally distributed random variable with mean 0 and standard deviation 10 * ωR₀. The term 10 * ωR₀ * η_nn[1] can be interpreted as a multiplicative deviation from the typical value tvR₀, and η_nn[1] is simply that deviation standardized. Using η_nn[1] both as the random effect on the initial value of Response and as an input to the MLP allows the model to learn how the initial value of Response affects the dynamics of the system.
The fix function is used to fix the random effects η_nn as non-dynamic inputs to the MLP, effectively creating an “individual” neural network iNN for each subject. Since the MLP has 5 inputs and η_nn is 3-dimensional, the output of fix will be another function with the remaining 2 inputs.
5 Fitting the model
Next, we fit the model to the training data. The fit function is used to estimate the parameters of the model using the training population trainpop. The initial parameter estimates are provided by init_params(model), and we use the MAP(FOCE()) method for fitting. The reason MAP is required is that under the hood, the regularization is implemented as a prior distribution on the parameters. Without MAP, the regularization would not be applied. FOCE() indicates that the first-order conditional estimation method will be used to approximate the marginal likelihood when fitting. We run the fit up to 150 iterations. Note that this step will take some time.
fpm = fit(
model,
trainpop,
init_params(model),
MAP(FOCE());
optim_options = (; iterations = 300),
)6 Evaluating the model on the training data
We can calculate the individual and population predictions using the fitted model on the training data using predict and visualize the predictions using plotgrid.
pred_traindata = predict(fpm; obstimes = 0:0.1:24);
plotgrid(pred_traindata)
After fitting, we can use inspect function to compute a number of useful statistics using the fitted model, including the empirical Bayes estimates of the random effects, the individual and population predictions, and the weighted residuals. The output of inspect can be passed to a number of plotting functions to visualize the results. The goodness_of_fit is a generally useful function that can be used to visualize multiple goodness of fit metrics.
ins = inspect(fpm)
goodness_of_fit(ins)[ Info: Calculating predictions. [ Info: Calculating weighted residuals. [ Info: Calculating empirical bayes. [ Info: Evaluating dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
Another standard diagnostic plot in pharmacometrics is the visual predictive check (VPC). The VPC compares the observed data with the simulated data from the fitted model. We can use the vpc function to compute the VPC metrics and then plot them using vpc_plot.
vpc_res = vpc(fpm, observations = [:Outcome], bandwidth = 3)vpc_plot(
vpc_res;
simquantile_medians = true,
axis = (xlabel = "Time (hours)", ylabel = "Outcome"),
)
7 Evaluating the model on the test data
A model can generally be considered successful if its predictions match the observations well for test data, not seen during the fit, and if the other goodness of fit diagnostics are satisfactory. We can predict for the test data using predict and visualize the predictions using plotgrid.
pred = predict(model, testpop, coef(fpm); obstimes = 0:0.1:24);
plotgrid(pred; ylabel = "Outcome (Test data)")
We can also plot the VPC for the test data using the vpc function. This will give us an idea of how well the model generalizes to unseen data.
vpc_res_test = vpc(model, testpop, coef(fpm); observations = [:Outcome], bandwidth = 3)vpc_plot(
vpc_res_test;
simquantile_medians = true,
axis = (xlabel = "Time (hours)", ylabel = "Outcome"),
)
8 New dosing regimen
Assuming the actual relationship between the drug and response is identifiable from the data, if we discovered this relationship, we should be able to use the previously fitted model to fit the data for new synthetic subjects under a dosing regimen that the model has never seen before in the training data. Let’s see what happens if we supply three low doses (0.3), wait a while and then supply two high doses (1.5).
To test this, we generate some synthetic test data using the synthetic_data function. The dosing regimen is specified using the DosageRegimen function.
dr2 = DosageRegimen(0.3, ii = 3, addl = 2)
dr3 = DosageRegimen(1.5, time = 25, ii = 8, addl = 1)
testpop2 = synthetic_data(
datamodel,
DosageRegimen(dr2, dr3),
p_data;
nsubj = 12,
obstimes = 0:2:48,
)We then use the fitted model to predict the response for the new subjects.
pred2 = predict(model, testpop2, coef(fpm); obstimes = 0:0.01:48);
plotgrid(pred2)
We can then overlay the data-generating model’s individual predictions of these test subjects to compare them against our previously fitted model. Note that no model fine-tuning is done here, we are simply using the fitted model to predict the response for the new subjects.
pred_truth = predict(datamodel, testpop2, p_data; obstimes = 0:0.01:48);
plotgrid!(
pred_truth;
pred = false,
ipred = (; color = Cycled(3), label = "DataModel ipred"),
)
We can see that the model has successfully learned some generalizable dynamics for the system, as evidenced by the DeepNLME model’s ability to predict the test data accurately, including dosing regimens that it was not trained on.
9 Conclusion
In this tutorial, we have demonstrated how to use DeepNLME to model data generated from an indirect response model. We have shown how to define a neural-embedded NLME model, fit it to the training data, plot fitting diagnostics, and then validate the model’s performance on test data. The model successfully learned a good dynamics model that was able to generalize to new dosing regimens.