Wolfram Blog
Robert Palmer

Drug Dose Selection Using Wolfram SystemModeler and Mathematica

October 25, 2012 — Robert Palmer, Applications Engineer

During the last decades, the development and use of therapeutic monoclonal antibodies (mAbs) have grown rapidly. Today, more than 30 different mAbs are successfully used in the clinic—playing important roles in treating complex diseases such as cancers and auto-immune disorders—and more than 200 are in clinical trials.

The history of mAbs has, however, not been without problems. In 2006, a first-in-human clinical trial of an mAb, aimed at treating leukemia and rheumatoid arthritis, went terribly wrong. Although the trial was run according to an approved protocol, all volunteers receiving the drug had severe inflammatory reactions and multiple organ failure. The tragic event shocked the medical community and highlighted a very important issue: how do you select a safe starting dose in first-in-human trials?

Now, as you may guess, the complete answer to this question is not an easy one. It’s also beyond the scope of this blog post. However, as a consequence of the dramatic happenings in 2006, the European Medicines Agency (EMEA) recently published new guidelines to address the issue of starting dose selection in first-in-human trials. Interestingly, the guidelines recommend that the use of modeling and simulation should play an integral part in the selection process, and in this post I thought we would study what such an approach might look like using Wolfram SystemModeler and Mathematica.

Drug Dose Selection Using Wolfram SystemModeler and Mathematica

The Minimum Anticipated Biological Effect Level

To make things a bit easier to follow, let’s set up an example.

First, to keep things focused, we’ll limit the example to mAbs. Due to their very specific nature and long half-lives, selecting starting doses for mAbs can be tricky, and tools that can aid this process are of added importance.

Second, let’s hypothesize that we have developed a new type of therapeutic mAb that we want to test in a clinical trial. The EMEA proposes that we should use a trial-starting dose that will result in the minimum anticipated biological effect level (MABEL). The MABEL is usually defined by 1) looking at the most sensitive biomarker associated with the drug’s mechanism of action, and 2) specifying how much this biomarker is allowed to change. For instance, let’s say that our mAb works by binding to and neutralizing a certain type of target molecule—a common mechanism of action of mAbs. The most sensitive biomarker is then typically the target occupancy (TO), that is, the relative number of target molecules bound by the mAb. To define our MABEL, we can say that we, for example, want no more than 10% TO in our first trial, and finding a recommended starting dose is then a question of finding a dose that would give us the 10% TO.

This is a question we should now try to answer using modeling and simulation.

Plot showing our minimum anticipated biological effect level (MABEL)

The Target Mediated Drug Disposition Model

OK, so let’s begin with the task at hand. The first thing we need to do is build a model that somehow relates the mAb dose to our biomarker. In pharmacology, such a model is usually referred to as a PK/PD model, where the PK means that the model describes the pharmacokinetics of the drug (i.e. the absorption, distribution, and elimination of the drug), and the PD means that the model also describes the pharmacodynamics of the drug (i.e. the relationship between the drug’s concentration and its effect on interesting biomarkers). In principle, numerous approaches exist for building PK/PD models, but in this example, we’ll make use of a target-mediated drug disposition (TMDD) model, which is a type of model frequently used to describe the PK/PD of mAbs.

By using components from the BioChem library, we can easily put together the TMDD model in SystemModeler and get something looking like this:

Model built in Wolfram SystemModeler

From the diagram, we can decipher that the model describes mAb PK by 1) a user-specified input describing the rate of appearance of the mAb in blood plasma, 2) a tissue compartment to describe non-specific mAb tissue binding or distribution, and 3) two possible mechanisms for mAb elimination: either by direct elimination or by target binding and subsequent degradation. Further, the PD of the mAb is described by the target synthesis, degradation, and target-mAb binding reactions.

Now, with the TMDD model, we could actually turn directly to the question at hand and search for the dose giving us our 10% TO. Before we do this, though, it’s always useful to run a test simulation, just to see how the model behaves. We can do this conveniently in Mathematica.

Getting to Know Our Model

The first thing we need to do to simulate our model is to define the input signal, that is, the mAb rate of appearance in plasma. The rate is dependent on two things: the total drug dose and the route of drug administration. Thus, let’s say that we want to use a dose of 50 mg mAb and administer it by constant intravenous infusion during 15 minutes:

mgDose = 50; minInfusionTime = 15;

The rate of appearance in plasma, during the infusion period, will then simply be:



And we can define our input signal as a Piecewise function having the value 10/3 up to 15 min and zero otherwise:

Defining our input signal as a Piecewise function

With the input defined, we can then run our simulation. Using WSMSimulate, we simulate the model for 20 days:


simTime = 20*60*24;

sim50mgDose = WSMSimulate["TMDDModel", {0, simTime},  WSMInputFunctions → {"rateOfAppearance" → inputSignal}];

We then extract some variables we want to study. Let’s look at the mAb and target concentrations in plasma:

mAbConc = sim50mgDose[{"plasmaCompartment.mAb.c"}, t];

targetConc = sim50mgDose[{"plasmaCompartment.target.c"}, t];

Finally we plot the result using a few plot options to tweak the appearance:

Plotting the result

Plots of mAb concentration in plasma and target concentration in plasma

These plots give us a first idea of what should happen when the drug is injected into the body. Looking at the upper plot, we see that the mAb plasma concentration should reach its peak instantly and then decline in a biphasic pattern as a result of mAb distribution and elimination. The second plot shows a reverse behavior: the concentration of free target molecules should quickly decline due to mAb-target binding, and then slowly return to baseline as the mAb concentration decreases.

With respect to our question, the lowest observed target concentration seen in the second plot is particularly interesting to look at, since this concentration is directly related to the maximum target occupancy. As we can see, the target concentration reaches almost zero at the start of the simulation, indicating that a very high percentage of the target molecules has been bound by antibodies.

To get a better picture of the actual TO, we can calculate it using the following function:

targetOccupancy[   targetConc_] := (1 - targetConc/baselineTargetConc)*100

baselineTargetConc = 4.5;

Apply the function to the target concentration:

to50mgDose = targetOccupancy[targetConc];

And then look at the result:

Plot[to50mgDose, {t, 0, simTime}, Evaluate@to50mgPlotOptions]

Target occupancy

Note that the maximum target occupancy is almost 100%, and consequently much higher than the 10% we aim for.

Finding Our Dose

OK, so now we know how the model behaves, and we also know that a 50 mg dose is definitely too high. Let’s move on and try to find a more correct dose, still using Mathematica.

Finding the dose corresponding to 10% TO could be done in several ways. For instance, we could use a trial and error approach, or we could create an objective function and use the Mathematica function FindMinimum. A more illustrative way, though, is to create a dose-response plot of the maximum TO.

To do this, we’ll start by defining a range of different doses to investigate. We use a logarithmic scale to cover both high and low values:

mgDoses = Table[10^x, {x, -2, 2, 0.1}]

We then slightly change the model input so it works with any dose, and set up a ParallelTable to simulate our different doses:

inputFunction[mgDose_] :=    Function[t,     Piecewise[{{mgDose/minInfusionTime, t < minInfusionTime}}, 0]];


simDoses =    ParallelTable[    WSMSimulate["TMDDModel", {0, simTime}, WSMInputFunctions → {"rateOfAppearance" → inputFunction[dose]}], {dose, mgDoses}];

Using the same procedure as when plotting the simulations of the 50 mg dose, we get a result looking like this:

targetConcDoses =   Map[#[{"plasmaCompartment.target.c"}, t] &, simDoses];

toDoses = targetOccupancy[targetConcDoses];

Plot[toDoses, {t, 0, simTime}, Evaluate@toDosesPlotOptions]

Target occupancy

As we can see, we have covered the complete range of different TOs.

Now, based on the result above, we need to identify the maximum TO in each simulation. The easiest way to do this is using FindMaximum:

Identifying the maximum TO in each simulation

We can finally plot the maximum TOs versus our doses:

plotData = Transpose[{mgDoses, maxTO}];

ListLogLinearPlot[plotData, maxTODosesPlotOptions]

Target occupancy

With this plot, we get a good overview of how the TO changes with different doses. Looking at the curve, we can quickly see that the dose corresponding to 10% TO should be approximately 0.5 mg. By interpolating the data, we can get a more exact answer:

iData = Transpose[{maxTO, mgDoses}];



According to our model, we should consequently use a dose of 0.46 mg as the starting dose in our trial. Of course, the result should only be seen as a prediction and, in case of a real scenario, be regarded in the light of other information collected during the pre-clinical phase of the drug development. Nevertheless, the information we get out of the model is valuable, and with a model in place we can easily go ahead and investigate more questions. What happens if we change the drug administration route? What happens if we design the antibody to bind more efficiently to its target molecule? Or how frequently do we need to administer the mAb to guarantee a certain minimum TO?

I will not go through the answers to these questions here, but feel free to download the model below and try to answer them yourself. A notebook with all the Mathematica code used in this blog post is included as well. The notebook also includes a sensitivity analysis that shows how sensitive the dose prediction is to changes in model parameters.

Download this ZIP file that includes the post as a CDF file, the model, and the code used above.

Posted in: SystemModeler
Leave a Comment

No Comments

Leave a comment



Translate this page: