Skip to content

TBMT42 Lab 3: NLME and PKPD Modelling

This is the third computer exercise in the course TMBT42 - Systems biology, Digital Twins and AI.

In this computer exercise you will apply non-linear mixed effect modelling (NLME), allometric based interspecies translations, and pharmacokinetic/pharmacodynamic (PKPD) modelling, to produce human drug treatment predictions from animal pre-clinical data. Human drug treatment predictions without human drug data are an important tool used in modern drug development and is a practical example of what mathematical modelling can be useful for. In this exercise you are the data analyst that hold the fate of a substrate under drug development, called Substrate X1 (SubX1). You are asked to analyse pre-clinical animal data and determine its potential for human applications. Everyone in your team is waiting for your analysis - Should we continue investing resources into the substrate, or is it time to look for other candidates?

Throughout this lab you will encounter different admonitions, which as a type of text box. A short legend of the different kinds of admonitions used can be found in the box below:

Different admonitions used (click me)

Background and Introduction

Useful information

Here you are tasked to do something

Reflective questions

This computer exercise is developed in MATLAB, however, the method can be implemented in any environment of your choosing. Although, please note that the practical instructions below will focus on implementation in MATLAB and the supplied code is written in MATLAB. If you choose to proceed with the computer exercise in another environment/software than MATLAB you will have to compile the tools and re-write the pre-made scripts yourself.

This computer exercises contains solutions and toolboxes that requires the following to be installed before you start: :

  • MATLAB running on your computer (or a computer in the computer hall)
  • A C-compiler installed in MATLAB

  • Downloaded and installed IQM tools in MATLAB

  • Suitable toolbox packages installed. You will need tools for optimisation algorithms, statistical testing, and random sampling. In MATLAB you'll want the following toolboxes installed:

    • "Global optimization Toolbox"
    • "Statistics and Machine Learning Toolbox"

All files needed in this example are located in the folder NLME_PKPD inside the Lab3 files. Copy these files to your MATLAB directory.

How to pass this computer exercise

To pass this computer exercise you need to follow the instructions herein, write a report, and upload the report to the TBMT42 course page on Lisam.

The report has to firstly include a short abstract/summary (max 400 words) which describes what you have done in this exercise and what your recommendation is for the future of the substrate that you are analysing. After the short abstract/summary you have to answers all the questions that are listed at the bottom of this page. The questions can either be answered as one main text or in a list format with both question and answer written. In your report, please include relevant figures that you have produced throughout this exercise.

Introduction

Decreasing levels of HormoneX has been linked to several health benefits. However, all previous drug candidates that reduces HormoneX that your team has produced has been rejected because it does not fulfil the required substrate-specific criteria to ensure that the potential future drug would be useful, safe, and affordable. Today you got sent an email from your teams’ project leader that reads:

Email with instructions for SubX1 analysis

" Hi,

Our latest SubX1 candidate seems promising! Could you please "run the numbers" to see if this is the candidate we should be moving forward with? In this email you will find the material for your analysis, see the attached NLME_PKPD folder.

If we are going to proceed with SubX1, we need an oral administrated human treatment that fulfill the following criteria:

  1. The toxicology team concluded that SubX1 plasma concentrations can never be higher than 600 umol/L.

  2. The biology and metabolism team concluded that HormoneX has to be reduced by 60% from the patient's basal levels for the health benefits to take effect.

  3. The production team concluded that we cannot move forward if the treatment requires more then 0,055 mol/week/patient.

  4. The marketing team concluded that the consumer is not interested in a pill that has to be taken more than twice a day.

We have a meeting with the board in a week please have your analysis ready for then.

Good luck,

project leader"

In the NLME_PKPD folder you will find:

  • "Data located in the folder Data, which includes a SubX1 rat experiment (N=10), SubX1 monkey experiment (N=1), and a SubX1 dog experiment (N=1). All experiments are in the same SubX1_data file. The data contains information on SubX1 and HormoneX plasma concentrations after an oral dose of 80 umol/kg in all species. In the data folder there is also an Excel file with the same information as the SubX1_data file."

  • "A mathematical model called SubX1_model.txt can be found in the main folder. The mathematical model was created by a colleague during the development of a previous HormoneX decreasing candidate.

  • The toolbox used to simulate the SubX1_model.txt is in the folder IQMtools which should be installed on your computer. If not installed, please read the installation.txt found in the IQMtools folder.

  • "Code as a foundation for your analysis is supplied according to a typical three-step method that your team is using: (i) Step1_PopulationVariability, (ii) Step2_MetabolicTranslation, and (iii) "Step3_ClinicalTesting.

There are many different ways you can solve your task. In summary, this exercise will exemplify the following solution:

  • Step 1: Find participant-specific variability from the rat data (N=10).

  • Step 2: Find metabolic scaling variables which describes a correlation between reaction rates and volumes to total species body weight. This is done using the rat, monkey, and dog experiments.

  • Step 3: Make human drug treatment prediction from the results of step 1 and 2.

In this three-step method there are many simplifications and assumptions which is interesting to discuss, explore, and/or improve. Simplifications and assumptions are a necessary part of the analysis and this exercise. You are encouraged to change these assumptions to your own hypothesis - if something doesn't work then you can always return to the original pre-written method.

Best of luck and have fun!

Getting familiar with the data and model

Open the main.m script. The main script runs all other scripts to ensure that they are done in the correct order. When running a script from the main it is important to remove the # symbol. Make sure to add the # symbol before a script in the main when you don’t want to run it. The first script that the main.m calls for is LoadThings.m, which is a script designed to always run when the main.m is used.

Open the LoadThings.m script and get familiar with what the script does. Everything you write here will be run every time you call on something through the main.m script. Typical things that should be included here is setting up your workspace, setting up the model you work with, and declaring variables that all scripts require, e.g. colours, index names, species body weight, etc.

As the first step you will be asked to plot and get familiar with the experiment data found in the Data folder (see Task 1 below). To achive this task you have been supplied with two scripts called PlotRatData.m and PlotMetabolicScalingData.m. These two scripts are located in the Step1_PopulationVariability and Step2_MetabolicTranslation folders. Both scripts are a foundation for Task 1 and you can choose if you want to use the scripts or make your own plots. Don't forget to remove the # symbol in the main if you want to call for the pre-made script.

Task 1: Data processing

Task 1A: Plot the pharmacokinetics of SubX1

The pharmacokinetic data of SubX1 is provided as a concentration in plasma (umol/L). The plots should include the rat experiment (N=10), monkey experiment (N=1), and dog experiment (N=1). One typical pharmacokinetic plot is to have SubX1 plasma concentration on the y-axis and time on the x-axis.

What is pharmacokinetics?

Pharmacokinetics includes drug rate of appearance, rate of disappearance, and distribution. It can be described/simplified as "what the system does to the drug".

Task 1B: Plot the pharmacodynamics of SubX1

The pharmacodynamic data of SubX is provided as HormoneX concentration in plasma (umol/L). The basal HormoneX concentration can be seen at 0 min. One possible solution to animals and species having different basal levels can be to normalise all subjects to their own baseline. To normalise metabolites loses information and can be a bad call when a saturation in reaction or transportation rates are expected. However, normalisation can also be a strong tool in some situations where there are differences in baselines and the aim is to reduce metabolites below patient-specific thresholds. The plots in this task should include the HormoneX levels from all three rat, monkey, and dog experiments. Your following analysis can be either be based on a normalised HormoneX or the concentration of HormoneX. The pre-prepared scripts and model is based on a normalised HormoneX levels.

What is pharmacodynamics?

Pharmacodynamics includes the effect of the drug. It can be described/simplified to "what the drug does to the system", which in this case is how HormoneX changes.

Task 1C: Visually compare your pharmacokinetic and pharmacodynamic plots

Look at the figures and consider the difference between the different rats and species. The dose is an oral administration of 80 umol/kg in all experiments, where the rats’ weight was [0.3,0.28,0.29,0.32,0.33,0.33,0.31,0.27,0.31,0.34] kg (mean: 0.308 kg, stdev: 0.023 kg), the monkey weight was 4.31 kg, and the dog weight was 9.75 kg. The weight of all subjects is included in the pre-written scripts for you. Please note differences in rate of drug clearance and drug distribution between species and rats. One effective tool can be to plot the data relative to species total body weight.

Reflective questions

  • Why do we see a rapid loss of plasma SubX1 in the first couple hours and then a lower loss of SubX1 closer towards the end of the experiment? Can you imagine how this would be mathematically described in your model?

  • Consider the correlation between SubX1 concentration and HormoneX plasma concentration. Can you see any difference in drug potency, drug clearance, and distribution?

After getting familiar with the data, it is time to get familiar with the SubX1_model.txt model that was sent to you:

States:

  • "State A is an "Administration state" which results in the delay between a drug SubX1 administration (umol) and when it is observed in the body. The administration delay can be used for many different administration routes, and in this example is used for oral administrations.

  • "State C1 A is "Central state 1" which represents drug concentration (umol/L) in the plasma"

  • "State C2 A is "Central state 2" which represents drug concentration (umol/L) not in the plasma, but still accounted for in the model. This includes substrate in organs and tissue."

  • "State HormoneX represents relative hormone levels to subject baseline (%) "

  • "State mod is a timer which can be used to re-administer drugs when the value of mod gets to the value ReAdministration_Time. This is declared in the MODEL EVENTS.

Parameteters:

  • Ka - rate of absorption (1/min)

  • Qref - rate of transportation between state C1 and C2 (1/min)

  • V1ref - volume of state C1 (L)

  • V2ref - volume of state C2 (L)

  • CLref - rate of drug loss (1/min)

  • ECmax - maximum drug effect relative to vehicle (ratio)

  • EC50 - drug potency constant (umol/L)

  • EChill - rate of effect slope around drug EC50 (dimensionless). The type of drug you are analysing has typically a value between 1 and 3. Setting a EChill value above 8 is an effective tool to describe switches and simultaneous effects.

  • Rb - Metabolic scaling constant for reactions (dimensionless)

  • Vb - Metabolic scaling constant for volumes (dimensionless)

  • BW - Total body weight (kg)

  • f - Bioavailability (%)

  • ReAdministration_Time - Time between re-administrations (min)

  • ReAdministration_Amount - Drug dose at re-administrations (umol)

How do SubX1 affect HormoneX in the model?

In the model, SubX1 affects HormoneX balance through increasing the negative derivative term of the state HormoneX. The increasing negative term represents a sum of all pathways a drug can be removed out of the subject, e.g. active clearance, excretion, degradation etc. The increased negative term is done through the MODEL REACTION SubX1_PD.

Parameter values and metabolic scaling The parameter values listed in the SubX1_model.txt model is used as a foundation to describe all species simulated in the model. Some reactions and all volumes are expected to be different between subjects and species. These differences can be described with random effects, which are steps from a population describing fixed effect, and/or fixed effect scaling to co-variables. The subject specific parameters (fixed effects + random effects) can be restricted to be normally distributed around the fixed effect when they are expected to be so in real life. Co-variables that scale fixed effects reduces the need of random effects and typically reduces model predictions uncertainties. Co-variables can be a Boolean, or variables such as total body weight, organ surface area, age, etc.. If the covariable is an allometric variable then the process is called Allometric scaling. A commonly used co-variable for the type of analysis you do in this task is an allometric scaling based on total body weight. The allometric scaling in the SubX1_model.txt model is done in the MODEL VARIABLES section.

Allometic scaling and Kleiber's law

Allometric scaling is standard practice in drug development research and clinical applications. An example of clinical application of allometric scaling in your everyday life is when your doctor prescribes a drug treatment dependent on your total body weight, this is typically done with a 1:1 ratio (dose/kg). In drug development research the allometric scaling used is typically around Ref(BWref/BWact)^1 for volumes and Ref(BWref/BWact)^0.75 for metabolic rates when the central compartment unit is concentration. In this scaling function, Ref is the parameter value describing the subject with the lowest bodyweight in the estimation data, BWref is the total bodyweight of the subject with the lowest bodyweight in the estimation data, and BWact is the bodyweight of the subject the function scales the parameter value to. The correlation starting guess of ^1 for volumes and ^0.75 for metabolic rates is called Kleiber's law.

Step 1

Now when you are familiar with the SubX1_model.txt model and the PlotRatData.m data, it is time to put them together and simulate the rat experiment with the model (Task 2). The shared model parameter values in the SubX1_model.txt describes a previous HormoneX reducing candidate and will not describe the new SubX1 candidate well.

Task 2: Model formulation and initial guess of rat experiment

Task 2A: Simulate the rat experiment with the model that you were supplied with.

To help with this task you have been supplied with the script PlotRatExperiment_InitialGuess which is found under the Step1_PopulationVariability folder. Don't forget to remove the # symbol in the main if you want to call for the pre-made script.

Reflective questions

  • Was the allometric scaling enough to describe all the differences seen in the rats?

  • In what scenario do you think the presence of allometric scaling based on total body weight is the most important? Remember the text above which explained how random effects and allometric scaling can be used together.

  • In what scenario do you think the presence of allometric scaling can make the model worse and increase model uncertainties?

Task 2b: Model structure and parameter initial guesses

Manually change model parameter in the SubX1_model.txt and see how it affects the model agreement to data. The parameter values listed in the model file will be the initial guess for when you are optimising the model to describe the data (Task 3). If there is anything you want to try and change in the model then from now on you are welcome to do so. Please keep in mind the important balance between model complexity and model agreement to data. Changing the model is completely optional but to test somethings and then potentially go back to the original model is encouraged to further increase your understanding of Systems biology. Please think twice before you make change and rename things, as the pre-existing scripts are written to work on the pre-supplied model.

Task 2C: Save the model initial agreement to the experiment

Save plots showing the agreement between model and the rat experiment with the unoptimized model parameters.

As seen in the rat experiment, there are fundamental and natural variability between the different subjects which is not fully described by their differences in the co-variable used in the allometric scaling functions (total body weight). Therefore, non-linear mixed effect modelling (NLME) can be used to add random effects to describe variability seen between subjects. The usage of a combined NLME and co-variables is a powerful tool and is a main point you should learn from this exercise.

Task 3: Train the model on the rat data and plot the agreement to data

Task 3A: Run a NLME optimisation to describe the rat experiment.

To succeed with this task you have been supplied with the script OptimiseRatExperiment which is found under the Step1_PopulationVariability folder. You are welcome to change the optimisation settings and NLME optimisation method if you want to. When the NLME optimisation finishes you will be provided with the outputs PHI, PSI, b, and PHIR which will be saved as a file named OptimisedParams_Step1. The variables PHI, PSI, b, and PHIR meaning can be derrived from the code combined with the nlmefit documentation found on MATLAB webpage and through writing help nlmefit in the MATLAB terminal. You can also look at the script PlotRatExperiment in the Step1_PopulationVariability folder to see how the different variables are used. Don't forget to remove the # symbol in the main if you want to call for the pre-made script.

Where is PHI, PSI, b, and PHIR created?

When the NLME optimisation script is finished then you will see PHI, PSI, b, and PHIR in your workspace. PHI, PSI, and b are outputs from the NLME optimisation and PHIR is calculated after the NLME optimisation from PHI and b.

Short description on how the OptimiseRatExperiment function works

The OptimiseRatExperiment calls for a NLME optimisation. The NLME optimisation loops over the SimulateRatExperiment function to simulate all the subjects and minimize the differences between the variable Y (found in the optimisation script) and the variable varargout (found in the simulation script). In the OptimiseRatExperiment the variable GROUP is declared, which is a vector that supplies the NLME optimisation with information of which datapoints is from which subject. In the shared NLME optimisation all subjects are supplied with unique random effects for the parameters which is allowed to have random effects. The vector REParamsSelect determines which parameters are allowed to have a random effect. All parameters not declared in the REParamsSelect vector will be restricted to only having one fixed effect value to describe the full population. The global variable GroupID is one solution to add allometric scaling to the NLME optimisaiton, where eatch rat is simulated with different total body weight (parameter BW). Please feel free to give more or less parameters variability and remember to add the variability to the fixed effects after the optimisation, see the calculation of PHIR to do so.

Understanding the NLME optimisation output
  • What does PHI, PSI, b, and PHIR represent? Please explain what they are and what they can be used for.
Task 3b: Plot the model agreement to data, parameter variability, and a PKPD curve after the NLME optimisation

To succeed with your task you have been supplied with the script PlotRatExperiment which is found under the Step1_PopulationVariability folder. Take a minute to compare the pharmacokinetics to the observed pharmacodynamics. Please specify what parameters have random effects and comment on the variability of the parameters.

PKPD curve

PKPD curve is a powerful tool which will help you to visualize drug potency and effect. The typical PKPD curve has pharmacodynamics the y-axis and drug concentration on the x-axis.

With a model describing the rat experiment (N=10) you have successfully found parameter variability which describes the population variability seen within the same species.

Step 2

You can now use the result from Step 1 to make predictions of species with other bodyweights. However, the allometric scaling will not be optimal without a calibration to a population with a larger variability in body weights. Before an allometric scaling calibration, the typical metabolic rates and volumes can be assumed to follow "Kleiber's law". To do this calibration your team founded a small monkey (N=1) and dog experiment (N=1), which can be combined with the rat data you worked with in Step 1.

With a calibrated substrate-specific allometric scaling the interspecies predictions will scale volumes and metabolic rates based on total body weight with a lower uncertainty than previously. However, it will not take into consideration other fundamental differences between different species. One of many fundamental differences between species is bioavailability.

Bioavailability and more: Not the same in all species

Bioavailability can be simplified and reduced to a percentage of how much of a drug that enters the body reaches the plasma. Bioavailability has no found correlation to total body weight, and thus it is not obvious on how you should take bioavailability into consideration in your analysis. In the shared scripts the bioavailability is assumed to be 62% for rats, 77% for monkeys, and 52% for dogs, which is the bioavailability from a previous HormoneX decreasing candidate with a similar molecular structure to SubX1. Typically, before human drug testing the human drug bioavailability is estimated from in vitro experiments and/or assumed to have a specific value based on other species bioavailability. This assumption could for example be a mean bioavailability of all species or the bioavailability of a species with the lowest difference to human digestive system. In the supplied code human bioavailability is assumed to be identical to the monkey bioavailability of the previous candidate (77%). One way of find the oral bioavailability of a drug is to compare IV administrations to an oral administration, where the IV administration by definition has a 100% bioavailability. However, IV administrated SubX1 data is not available to you, and you will be forced to make an assumption. Similar assumptions can be made with drug potency, and metabolic rates when needed.

Don't forget to add the # symbol in the main on the scripts you don’t want to use. Re-running the optimisation in step 1 could yield new results which will overwrite your old Step 1 optimisation results.

Task 4: Determining the metabolic variable scaling constants

Task 4A: Predict the monkey and dog data based on your results from Step 1.

To succeed with your task you have been supplied with the script PlotMetabolicScaling_InitialGuess which is found under the Step2_MetabolicTranslation folder. The PlotMetabolicScaling_InitialGuess will help you plotting an initial guess of the monkey and dog experiment. Comment on the agreement to data and try to visualise/guess the variability needed to describe the monkey and dog data.

Task 4b: Optimise the model while taking into consideration all rat, monkey, and dog data

To suceed with your task you have been supplied with the script OptimiseMetabolicScaling which is found under the Step2_MetabolicTranslation folder. The optimisation in this task is similar to what was done in the Step 1 and please feel free to change the scripts and functions to your own liking. Changing the optimisation method is completely optional, however, it is encouraged to test something and then potentially go back to the original method if needed.

Task 4C: Plot the model agreement to data, plot the PKPD curve, and analyse the parameter variability.

To succeed with your task you have been supplied with the script PlotMetabolicScaling which is found under the Step2_MetabolicTranslation folder. Compare the results from this optimisation and the optimisation in Step 1.

Reflective questions

  • Is the parameter variability similar to the results seen in step 1?

  • What are the allometric scaling variable values Rb and Vb? Remember what they represent! Why do we have two of both of them?

  • Do the PKPD curve look as you expected and is the difference in drug effect between the subjects only because of differences in pharmacokinetics?

You have now used the rat data to quantify natural population-based variability (Step 1) and found metabolic reaction rate correlation to total species body weight (Step 2), you are now ready for the final step which is making a human drug prediction without human drug data (Step 3).

Step 3

You are now close to your end analysis and your project team is waiting anxiously for your conclusion. In the folder Step3_ClinicalTesting you will find the script PredictClinicalTesting which will help you in making the human drug treatment prediction. Look over the PredictClinicalTesting script, understand what it does, and make any potential changes to the prediction method. In the script feel free to changing the HumanBW variable to the total body weight you want to simulate and play around with the ReAdministration_Time variable to choose the drug administration frequency (minutes in the original model), and the ReAdministration_Amount for the amount of active substrate (SubX1) in the drug (umol).

Task 5: Test different treatments and save an acceptable treatment or the best treatment you can find.

Save a treatment that fulfils all the criteria listed in the instructions (see introduction section above). If your prediction of SubX1 in human cannot fulfil the criteria then save the treatment that you can produce that is the closest to the criteria. Do you believe SubX1 is good enough to invest more resources into? Your analysis quality is measured in the plausibility of your assumptions and the motivation you make regarding your predictions. How well you succeed as a pharmaceutical analytics is independent if the candidate is rejected or not, but will be affected if you reject a good candidate or recommend to continue with a bad candidate. Don't forget to remove the # symbol in the main if you want to call for the pre-made script.

What are the assumptions made to make the human prediction?

  • What assumptions have you made to produce the human predictions?

  • How did you use the rat, monkey, and dog data to produce the human prediction?

  • What human drug potency (EC50) is used and what is it derived from?

  • What human drug bioavailability (f) is used and what is it derived from?

  • What human drug drug clearance (CL) is used and what is it derived from?

  • What human drug plasma volume (V1) is used and what is it derived from?

You are now done with the modelling work! Well done!

Your final task is to write your analysis before the big board meeting (instructions below).

What should your computer exercise report include?

Your report should include a short abstract/summary and answer the questions below.

The abstract/summary should include a summary of what you have done and what your recommendation is for the future of SubX1 (max 400 words).

After the abstract/summary a text should follow where you answer all the questions below either as a list (question and then answer), or as a full main text. Make sure to include relevant figures you have produced during this exercise as a foundation for your analysis:

Task 1 - Data processing

  • What is pharmacodynamics and pharmacokinetics?

  • Why do you we see a quicker loss of plasma SubX1 in the first couple hours compared to the later datapoints?

  • How would you simulate an IV administration of SubX1 in your model?

Task 2 - Model formulation and initial guess of rat experiment

  • Was the allometric scaling enough to describe all the differences seen in the rats?

  • If you can only choose to use either allometric scaling or random effects, when would it be best to use the allometric scaling and when would it be best to only use random effects?

Task 3 - Train the model on the rat data and plot the agreement to data

  • What is the PHI, PSI, b, and PHIR variable after the NLME optimisation?

  • Shortly describe the optimisation method you used.

  • Show the model agreement to data after the Step 1 optimisation, parameter values, parameter variability, parameter correlation, and the PKPD curve.

Task 4 - Determining the metabolic variable scaling constants

  • Did the model agreement to the rat data change when you included the monkey and dog data?

  • What was the allometric scaling variable values Rb and Vb values you got from the optimisation? If you don’t have enough information to identify the Rb and Vb values, what is a typical initial guess you can use? What are Rb and Vb used for?

  • In your work, which parameters were given random effects, which parameters were scaled with allometric scaling, and which were described using both?

  • Show the model agreement to data after the Step 2 optimisation, parameter values, parameter variability, parameter correlation, and the PKPD curve.

Task 5 - Produce the final human treatment prediction of SubX1

  • What is your final decision on the future of SubX1? Describe the treatment you produced that either fulfilled the criteria (see introduction) or was the closes to fulfilling the criteria?

  • Did you find a significant difference in taking the same amount of substance either once a day or separated into smaller administrations? Did you expect those results?

  • List a few assumptions you made to produce the human treatment prediction. For example, what parameter values did you use for the prediction (EC50, f, CL, and V1)?

  • Are there any future experiments for SubX1 that you want to recommend before human testing?