NLMEModeling - 
Nonlinear Mixed Effects Modeling of Dynamical Systems

Mats Jirstrand, Fraunhofer-Chalmers Centre

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

Nonlinear mixed effects (NLME) modeling is a powerful tool to analyze time-series data from several individual entities. In this talk, we will give a brief overview of a package for NLME modeling in Mathematica entitled NLMEModeling, implementing the so called first-order conditional estimation method with sensitivity equation-based gradients for parameter estimation. NLMEModeling supports mixed effects modeling of dynamical systems where the underlying dynamics are described by either ordinary or stochastic differential equations combined with observation expressions with flexible observation error models. Moreover, NLMEModeling is a user-friendly package with functionality for parameter estimation, model validation (such as goodness-of-fit analysis and visual predictive checks), and generation of synthetic data by model simulation.

The package is freely available and provides an extensible add-on for Mathematica users facing the problem of estimating dynamical models  given time-series data from multiple individuals sharing the same model structure but with parameters from a population distribution. The package has been applied to challenging problems both in pharmacometrics (clinical and pre-clinical data) as well as in microbiology (multiple single cell data).

What is Nonlinear Mixed Effects Modeling?

Nonlinear Mixed Effects Modeling is a parameter estimation methodology very useful when model parameters vary between individual subjects and amount of data per subject are sparse

Nonlinear = nonlinear in the parameters

Mixed effects = fixed and random effects (‘effects = parameters’)

Fixed effects = parameters common to all subjects

Random effects = parameters assumed to vary between individuals according to a probability distribution

What is Nonlinear Mixed Effects Modeling?

What is Nonlinear Mixed Effects Modeling?

A Pharmacokinetic Model (ODE)

NLME Modeling Data

Let’s consider a one-compartment pharmacokinetic system with exponentially decreasing absorption

The system description (including initial conditions and an equation for an auxiliary variable)

sys = {
A1'[t]==-ka*A1[t], A1[0]==Dose,
A2'[t]==ka*A1[t]-CL/V*A2[t], A2[0]==0,
c[t]==A2[t]/V
};

Let’s have a look at the solution

sol=ParametricNDSolveValue[sys,c[t],{t,0,24},{ka,CL,V,Dose}]

for a given set of parameter values ((Subscript[k, a])=1, CL=10, V=50) and Dose=100.

Plot[sol[1,10,50,100],{t,0,24}]

Now, introduce some “random effects” (Subscript[η, 1])~N(0,0.3) to describe the distribution of parameters in a population of individuals

sys = {
A1'[t]==-ka*A1[t], A1[0]==Dose,
A2'[t]==ka*A1[t]-CL*Exp[η1]/V*A2[t], A2[0]==0,
c[t]==A2[t]/V
};

compute the solution

sol=ParametricNDSolveValue[sys,c[t],{t,0,24},{ka,CL,V,Dose,η1}]

sample 20 individuals and visualize their pharmacokinetic responses

SeedRandom[123];
plPK=Plot[Evaluate[Table[sol[1,10,50,100,η1],{η1,η1pop=RandomVariate[NormalDistribution[0,0.3],20]}]],{t,0,24},PlotRange->{-.1,2.1}]

We may also measure each individual according to some pre-specified sampling protocol, let’s say

tSampleTimes={0.25,0.5,1,1.5,2,3,4,6,8,12,18,24};

and assume a normally distributed measurement error proportional to the measured concentration: y((Subscript[t, k]))=c((Subscript[t, k]))+c((Subscript[t, k])) e((Subscript[t, k])), e((Subscript[t, k]))~N(0,0.2).

data=Table[{t,#+#*RandomVariate[NormalDistribution[0,0.2]]},{t,tSampleTimes}]&/@Table[sol[1,10,50,100,η1],{η1,η1pop}];

ListLinePlot[data,PlotMarkers->Automatic]

Using the NLMEModeling application

Load the package

Needs["NLMEModeling`"]

System equations and measurement model

sys = {
A1'[t]==-ka*A1[t], A1[0]==Dose,
A2'[t]==ka*A1[t]-CL*Exp[η1]/V*A2[t], A2[0]==0,
c[t]==A2[t]/V
};
obs={c[t]};

Covariance of the measurement noise and population model

realΣ={{(0.2*c[t])^2}};
realΩ={{0.3^2}};

Experiment with 20 individuals.
Define experimental protocol (sample times, dose). 
Set up data structure to generate synthetic population time-series data.

SeedRandom[31];
nInd=20;
timePoints={0.25,0.5,1,1.5,2,3,4,6,8,12,18,24};
dataStructure=Table[{timePoints,{SubjectID->ToString[i],Dose->100,
Weight->Round[RandomVariate[NormalDistribution[70,5]]]}},{i,1,nInd}]

Generate population data

{randomParamValuesData,data}=NLMEDynamicalModelSimulate[dataStructure,{sys,obs},
{{ka,1},{CL,10},{V,50}},{η1},Sigma->realΣ,Omega->realΩ]

and visualize the result

dataPlots=PlotNLMEData[data]

pkPlot=Show[dataPlots[[1,1]],FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic},Frame->{True,True,False,False},FrameLabel->{"Time after dose (h)","Drug concentration (mg/L)"}]

Estimate model utilizing the synthetic data

modelFit=NLMEDynamicalModelFit[data,{sys,obs},{{ka,0.9},{CL,5},{V,30}},η1,
Sigma->"Proportional"]

Export fitted model object for later use ...

(* Export[NotebookDirectory[]<>"Example1-FittedNLMEModel.m",modelFit] *)

Model Validation

Import a fitted model object

modelFit = Import[NotebookDirectory[]<>"Example1-FittedNLMEModel.m"]

The properties of the FittedNLMEModel object

modelFit["Properties"]

A model summary

modelFit["ModelSummary"]

We can plot the population model fit

PlotPopulationModelFit[modelFit][[1,1,1]]

and the individual model fits (i=1,2,3)

PlotIndividualModelFits[modelFit][[1,1;;3]]

Basic goodness-of-fit plots

GoodnessOfFitAnalysis[modelFit][[1,1]]

Visual predictive check (VPC) is a graphical visualization comparing empirical percentiles of the observed data with percentiles of simulated data obtained by Monte Carlo simulations of the model to be evaluated and the experimental study design behind the observed data.

plVPC=VisualPredictiveCheck[modelFit,200,Quantiles -> {0.1,0.5,0.9},ConfidenceInterval->90];

Show[plVPC,FrameLabel->{"Time after dose (h)","Concentration (mg/L)"},FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic},Frame->{True,True,False,False}]

Pharmacokinetic-Pharmacodynamic Model (ODE)

Let’s consider a one-compartment pharmacokinetic system with exponentially decreasing absorption, which drives a pharmacodynamic response modeled using a turn over model with input saturation.

System specification and observation expressions

sys = {
A1'[t]==-ka*A1[t],
A2'[t]==ka*A1[t]-phi1/V*A2[t],
R'[t]==kout*(phi2*(1-1*c[t]/(EC50+c[t]))-R[t]),
A1[0]==Dose,
A2[0]==0,
R[0]==phi2,
phi1==CL*Exp[η1],
phi2==E0*Exp[η2],
c[t]==A2[t]/V
};
obs={c[t],R[t]};

Synthetic Data Generation

Load the package

Needs["NLMEModeling`"]

Observation covariance matrix and population covariance

realΣ={{0.1^2+(0.2*c[t])^2,0},{0,5^2}} ;
realΩ={{0.3^2,0},{0,0.1^2}};

Define number of subjects (45), sampling times and covariates (e.g., dose 100, 300, 1000)

SeedRandom[100];
nGroup=15;
timePoints={0.25,0.5,1,1.5,2,3,4,6,8,12,18,24};
doseGroup1=Table[{timePoints,{SubjectID->ToString[i],Dose->100,
Weight->Round[RandomVariate[NormalDistribution[70,5]]]}},{i,1,nGroup}];
doseGroup2=Table[{timePoints,{SubjectID->ToString[i],Dose->300,
Weight->Round[RandomVariate[NormalDistribution[70,5]]]}},{i,nGroup+1,2*nGroup}];
doseGroup3=Table[{timePoints,{SubjectID->ToString[i],Dose->1000,
Weight->Round[RandomVariate[NormalDistribution[70,5]]]}},{i,2*nGroup+1,3*nGroup}];
dataStructure = Join[doseGroup1,doseGroup2,doseGroup3]

Generate dataset

SeedRandom[43];
{randomParamValuesData,data}=NLMEDynamicalModelSimulate[dataStructure,{sys,obs},
{{ka,1},{CL,10},{V,50},{kout,0.4},{E0,100},{EC50,2}},
{η1,η2},Sigma->realΣ,Omega->realΩ];

Export dataset

(* Export[NotebookDirectory[]<>"Example2-data.m",data]; *)

Data Exploration

Import data

data=Import[NotebookDirectory[]<>"Example2-data.m"];
dataPlots=PlotNLMEData[data];

Visualize PK

pkPlot=Show[dataPlots[[1,1]],Frame->{True,True,False,False},FrameLabel->{"Time after dose (h)","Drug concentration (mg/L)"},FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic}]

Visualize PD

pdPlot=Show[dataPlots[[1,2]],Frame->{True,True,False,False},FrameLabel->{"Time after dose (h)","Drug response, R(t)"},FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic},PlotRange->{Automatic,{0,140}}]

Parameter Estimation

Load the package

Needs["NLMEModeling`"]

Load data to be used

data=Import[NotebookDirectory[]<>"Example2-data.m"];

System equations and measurement expressions

sys = {
A1'[t]==-ka*A1[t],
A2'[t]==ka*A1[t]-phi1/V*A2[t],
R'[t]==kout*(phi2*(1-1*c[t]/(EC50+c[t]))-R[t]),
A1[0]==Dose,
A2[0]==0,
R[0]==phi2,
phi1==CL*Exp[η1],
phi2==E0*Exp[η2],
c[t]==A2[t]/V
};
obs={c[t],R[t]};

Measurement noise covariance matrix

SigmaMatrix=(({{add1^2+(prop1*c[t])^2, 0}, {0, add2^2}}));

Estimate parameters

modelFit=NLMEDynamicalModelFit[data,{sys,obs},
{{ka,0.8},{CL,20},{V,40},{kout,0.2},{E0,90},{EC50,5}},{η1,η2},
Sigma->{"Advanced",{{SigmaMatrix},{{add1,0.1},{prop1,0.1},{add2,6}}}},
Omega->"Full"]

Export fitted model object for later use ...

(* Export[NotebookDirectory[]<>"Example2-FittedNLMEModel.m",modelFit] *)

Model Validation

Import the fitted model object

modelFit=Import[NotebookDirectory[]<>"Example2-FittedNLMEModel.m"]

Properties of FittedNLMEModel-objects

modelFit["Properties"]

A model summary

modelFit["ModelSummary"]

We can plot the population model fit

(* PlotPopulationModelFit[modelFit] *)

and some individual model fits (i=1,2,3)

PlotIndividualModelFits[modelFit][[1,1;;3]]//GraphicsGrid

Basic Goodness-of-fit plot

{{pkGof,pdGof}}=GoodnessOfFitAnalysis[modelFit];

Pharmacokinetics (drug concentration) GOF-plot

pkGof

Pharmacodynamics (drug effect) GOF-plot

pdGof

Analysis of empirical Bayes estimates of the random effects (qq-plots, pair-plots, correlations)

EBECorrelationAnalysis[modelFit]

Box and Whisker charts for the empirical Bayes estimates of the random effects

EBEBoxWhiskerChart[modelFit]

We may also check for patterns in the EBEs related to different covariates (here: stratifying the EBEs wrt dose)

EBEBoxWhiskerChart[modelFit,Stratify->"Dose"]

Visual predictive check

VPCPlots=VisualPredictiveCheck[modelFit,200,
Quantiles -> {0.1,0.5,0.9},
ConfidenceInterval->90,
PredictionCorrection-> True
];

{pkVPC,pdVPC}=VPCPlots;

VPC for PK-model

pkOut=Show[pkVPC,Frame->{True,True,False,False},FrameLabel->{"Time after dose (h)","Drug concentration (mg/L)"},FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic}]

VPC for PD-model

pdOut=Show[pdVPC,Frame->{True,True,False,False},FrameLabel->{"Time after dose (h)","Drug response, R(t)"},FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic},PlotRange->{Automatic,{0,140}}]

Pharmacokinetic-Pharmacodynamic Model (SDE)

Let’s consider a one-compartment pharmacokinetic system with exponentially decreasing absorption, which drives a pharmacodynamic response modeled using a turn over model with input saturation. In addition we have an unspecified system noise term on the response dynamics.

System specification and observation expressions

sys = {
A1'[t]==-ka*A1[t],
A2'[t]==ka*A1[t]-phi1/V*A2[t],
R'[t]==kout*(phi2*(1-1*c[t]/(EC50+c[t]))-R[t])+g*w[t],
A1[0]==Dose,
A2[0]==0,
R[0]==phi2,
phi1==CL*Exp[η1],
phi2==E0*Exp[η2],
c[t]==A2[t]/V
};
obs={c[t],R[t]};

Synthetic Data Generation

Load the package

Needs["NLMEModeling`"]

Observation covariance matrix and population covariance

realΣ={{0.1^2+(0.2*c[t])^2,0},{0,5^2}} ;
realΩ={{0.3^2,0},{0,0.1^2}};

Define number of subjects (45), sampling times and covariates

SeedRandom[101];
nGroup=15;
timePoints={0.25,0.5,1,1.5,2,3,4,6,8,12,18,24};
doseGroup1=Table[{timePoints,{SubjectID->ToString[i],Dose->100,
Weight->Round[RandomVariate[NormalDistribution[70,5]]]}},{i,1,nGroup}];
doseGroup2=Table[{timePoints,{SubjectID->ToString[i],Dose->300,
Weight->Round[RandomVariate[NormalDistribution[70,5]]]}},{i,nGroup+1,2*nGroup}];
doseGroup3=Table[{timePoints,{SubjectID->ToString[i],Dose->1000,
Weight->Round[RandomVariate[NormalDistribution[70,5]]]}},{i,2*nGroup+1,3*nGroup}];
dataStructure = Join[doseGroup1,doseGroup2,doseGroup3];

Generate dataset

SeedRandom[44];
{randomParamValuesData,data}=NLMEDynamicalModelSimulate[dataStructure,{sys,obs},
{{ka,1},{CL,10},{V,50},{kout,0.4},{E0,100},{EC50,2},{g,10}},
{η1,η2},w,Sigma->realΣ,Omega->realΩ];

Export dataset

(* Export[NotebookDirectory[]<>"Example3-data.m",data] *)

Data Exploration

Import data

data=Import[NotebookDirectory[]<>"Example3-data.m"];
dataPlots=PlotNLMEData[data];

Visualize PK

pkPlot=Show[dataPlots[[1,1]],Frame->{True,True,False,False},FrameLabel->{"Time after dose (h)","Drug concentration (mg/L)"},FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic}]

Visualize PD

pdPlot=Show[dataPlots[[1,2]],Frame->{True,True,False,False},FrameLabel->{"Time after dose (h)","Drug response, R(t)"},FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic},PlotRange->{Automatic,{0,140}}]

Parameter Estimation

Load the package

Needs["NLMEModeling`"]

Load data to be used

data=Import[NotebookDirectory[]<>"Example3-data.m"];

System equations and measurement expressions

sys = {
A1'[t]==-ka*A1[t],
A2'[t]==ka*A1[t]-phi1/V*A2[t],
R'[t]==kout*(phi2*(1-1*c[t]/(EC50+c[t]))-R[t])+g*w[t],
A1[0]==Dose,
A2[0]==0,
R[0]==phi2,
phi1==CL*Exp[η1],
phi2==E0*Exp[η2],
c[t]==A2[t]/V
};
obs={c[t],R[t]};

Measurement noise covariance matrix

SigmaMatrix=(({{add1^2+(prop1*c[t])^2, 0}, {0, add2^2}}));

Estimate parameters

modelFit=NLMEDynamicalModelFit[data,{sys,obs},
{{ka,0.8},{CL,20},{V,40},{kout,0.2},{E0,90},{EC50,5},{g,5}},{η1,η2},w,
Sigma->{"Advanced",{{SigmaMatrix},{{add1,0.1},{prop1,0.1},{add2,6}}}},
Omega->"Diagonal"]

Export fitted model object for later use ...

(* Export[NotebookDirectory[]<>"Example3-FittedNLMEModel.m",modelFit] *)

Model Validation

Import model object

modelFit=Import[NotebookDirectory[]<>"Example3-FittedNLMEModel.m"]

Properties of FittedNLMEModel-objects

modelFit["Properties"]

A model summary

modelFit["ModelSummary"]

Visual predictive check

SeedRandom[99];
VPCTab=Table[VisualPredictiveCheck[modelFit,200,
Quantiles -> {0.5},
ConfidenceInterval->90,
PredictionCorrection-> False,
Stratification->{True,{Dose->dosei}}
],{dosei,{100,300,1000}}];

VPC for PK-model

pkPlots=Table[Show[VPCTab[[i,1]],PlotRange->{Automatic,Automatic},AxesOrigin->{0,0},Frame->{True,True,False,False},FrameLabel->{"Time after dose (h)","Drug concentration (mg/L)"},FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic}],{i,1,3}]

VPC for PD-model

pdPlots=Table[Show[VPCTab[[i,2]],PlotRange->{Automatic,{0,140}},AxesOrigin->{0,0},Frame->{True,True,False,False},FrameLabel->{"Time after dose (h)","Drug response, R(t)"},FrameTicks->{{0,1,2,3,4,6,8,12,18,24},Automatic}],{i,1,3}]

NLMEModeling - Overview

Parameter Estimation

NLMEDynamicalModelFit - estimate a nonlinear mixed effects model

FittedNLMEModel - symbolic representation of a model

Properties: "ModelSummary", "FixedEffects", "OmegaMatrix", "SigmaMatrix", "BestFit", "ParameterCovarianceMatrix", "ObjectiveFunctionValue", ...

Model Validation

GoodnessOfFitAnalysis - goodness-of-fit plots for an estimated model

EBECorrelationAnalysis - correlation analysis of random effects empirical Bayes estimate

EBEBoxWhiskerChart - box-and-whisker chart of random effects empirical Bayes estimate

VisualPredictiveCheck - visual predictive check for an estimated model

Data Generation

NLMEDynamicalModelSimulate - generates population data for a model

Visualization

PlotNLMEData - plots population data

PlotIndividualModelFits - plots individual model predictions

PlotPopulationModelFit - plots population model prediction

PlotObservationsPredictions - plots observations vs population and individual predictions

PlotWeightedResiduals - plots weighted residuals vs individual predictions and time

Nonlinear Mixed Effects Modeling - Theory

Acknowledgements

Joachim Almquist, AstraZeneca (FCC alumni)

Helga Kristin Olafsdottir, Chalmers (FCC alumni)

Jacob Leander, AstraZeneca (FCC alumni)

Anna Johnning,  FCC

Julia Larsson, AstraZeneca (FCC alumni)

Felix Held, Chalmers (FCC alumni)

Simon Berglund Watanabe (FCC alumni)

Swedish Foundation for Strategic Research

See also

Leander, J. et al (2020), NLMEModeling: A Wolfram Mathematica Package for Nonlinear Mixed Effects Modeling of Dynamical Systems https://arxiv.org/pdf/2011.06879.pdf ()
Leander, J. et al (2021) Nonlinear Mixed Effects Modeling of Deterministic and Stochastic Dynamical Systems in Wolfram Mathematica https://doi.org/10.1016/j.ifacol.2021.08.394 ()
Jirstrand, M. (2021) https://www.youtube.com/watch?v=5SG1o0oCfRE ()
NLMEModeling – A Mathematica package for NLME modeling
http://www.fcc.chalmers.se/software/other-software/nlmemodeling/ ()