Pharmacokinetics and Pharmacodynamics

Overview

This app allows exploration of a model that is very similar to the Antiviral Treatment Model (which in turn is based on the Basic Virus Model.)

The new feature in this app is that the drug is explicitly modeled. Such models go by the name of pharmacokinetics/pharmacodynamics (PK/PD) models. Read about the model in The Model tab. Then, work through the tasks described in the What To Do tab.

Learning Objectives

The Model

Model Overview

This model consists of 4 compartments modeling the basic dynamics of a viral infection in the presence of a drug. In this model, we track the following entities by assigning each to a compartment:

In addition to specifying the compartments of a model, we need to specify the processes governing the model dynamics. For our system, we specify the following processes/flows:

  1. Uninfected cells are produced at some rate n and naturally die at some rate dU.
  2. Virus infects cells at rate b, with unit conversion factor g.
  3. Infected cells produce new virus at rate p and die at rate dI.
  4. Free virus is removed at rate dV or goes on to infect further uninfected cells.
  5. Once treatment has started (t > txstart), a new drug dose, C0, is administered at fixed time intervals, tinterval, which increases C at that time to C + C0. This is modeled by an instantaneous increase. In between every dose, drug kinetics is modeled as decaying exponentially at rate dC.
  6. The efficacy of the drug, e, depends on the concentration of the drug, given by a so called Emax-equation. Emax is the maximum efficacy at high drug concentration (and should be between 0 and 1), C50 is the drug concentration at which the drug has half its maximum efficacy, and the parameter k dictates how rapidly efficacy increases as drug concentration increases.

Model Diagram

The diagram illustrating this compartmental model is shown in the figure.

Model Diagram

Model Diagram

Model Equations

Implementing this model as a continuous-time, deterministic model leads to the following set of ordinary differential equations:

\[\begin{align} \dot U & = n - d_U U - bUV \\ \dot I & = bUV - d_I I \\ \dot V & = (1-e)pI - d_V V - gb UV \\ \dot C & = - d_C C, \qquad C=C+C_0 \textrm{ at } t = t_{interval} \end{align}\]

tinterval is the time at which a new drug dose is given. Prior to treatment start, i.e., t < txstart, there is no drug and C = 0. The drug efficacy e is given by its own equation which depends on C as follows:

\[e = E_{max} \frac{C^k}{C^k+C_{50}}\]

What To Do

The model is assumed to run in units of days.

Task 1

For the first few tasks, we consider an acute viral infection and treatment (e.g. influenza and neuraminidase inhibitor drugs). Set number of uninfected cells to 105, 10 virions, no infected cells. We make the assumption that on the timescale of an acute infection (several days), the processes of natural, uninfected cell turnover are so slow that they can be ignored. Set values for the uninfected cell birth and death rates to reflect this assumption. Assume that infected cells have an average lifespan of 1 day, virus of 12 hours. Set virus production rate to 10, infection rate to 10-5 and conversion factor to 1. Start simulation at time 0 and run for 20 days. For the drug, assume treatment starts at day 10 and occurs daily thereafter at a dose of 1. Assume the drug decays at a rate of 1 per day. Set the concentration-dependent drug effect parameter to 1, the half maximum effect level to 1 and max drug efficacy to 0. Run the simulation. You should get 28694 infected cells at the peak and a max drug dose of 1.57. Since the drug levels are fairly low, you might need to plot the y-axis on a log scale and use plotly and zoom into the drug part to see it’s up and down pattern.

Record

Task 2

Since we just had the maximum drug efficacy at 0, the drug had no impact at reducing virus production. Now, change to maximum impact by setting this parameter to 1. What do you expect to see? What do you see? You should notice only a small impact with the final number of uninfected cells having increased a bit.

Now, start drug treatment earlier, at days 8, 6, 4, 2, 0. Observe how the drug now has more of an impact. Not surprisingly, the earlier the drug is started, the more of an impact it has. The drug is however not strong enough to suppress the infection. We explore this further next.

Record

Task 3

Let’s explore the PK part of the model, i.e., the kinetics of the drug as described by the drug concentration equation. Reset all inputs. Then, set maximum drug efficacy to 1 and start drug on day 2. Set the plot engine to plotly with a log y-axis. Run the model.

By reading off the values for the drug from the graph (or looking at the reported values below the graph), convince yourself that the peak drug concentration starts at 1 and slowly over the first few days increases to around 1.5. Now set the drug decay rate to 5. Convince yourself that the peak drug concentration does not go beyond 1, i.e. there is no build-up. Now, set dC = 0.1. You’ll now see a build up to a concentration of almost 9.

Next, set drug decay rate back to dC = 1 and change treatment interval to 3 days. Run simulation and convince yourself that the drug concentration again does not really go above the administered dose of 1. Re-do for a decay rate of dC = 0.1. You’ll see some drug build-up, but less. You can keep exploring how dose, decay rate of drug and time between treatment impact the curve you see for C. The impact of each of those components should be fairly easy to understand, but it’s still worth playing with it a bit.

Record

Task 4

Now, we’ll explore the PD part of the model, i.e., the equation \(e = E_{max} \frac{C^k}{C^k+C_{50}}\), which describes the impact of the drug for a given drug concentration. Reset all inputs. Start treatment at day 2. Run the model for Emax = 0, you should get the same values as in task 1.

Now, explore Emax. Set it to 0.5, then 1. You should see a substantial increase in the impact of the drug, evidenced by higher uninfected number of cells remaining. You might habe noticed that by increasing Emax, the infection is delayed, which is especially noticable for Emax = 1. To account for this, set the simulation time to 100 days, then run again for the 3 different Emax levels (0, 0.5, 1). You’ll still find that an increased drug efficacy leads to fewer infected cells. Note that for the way we set up the model, the highest biologically reasonable value for Emax is 1, which means e = 1. While you can stick higher values into the model, they don’t make any sense (you’d get negative virus production). Depending on how the model is formulated, Emax values above 1 could be possible.

Next, explore the parameter k by running the simulation for values of 1, 2, and 3. Leave everything else as it is (i.e. Emax = 1, tfinal = 100). You’ll find that it doesn’t have much of an impact on the infection. Finally, we’ll explore the C50 parameter.

Set k back to 1, and run the model for C50 values of 1, 0.1 and 10. You’ll notice that this parameter has a strong impact.

Note that for all the changes you just did, the drug time-series (the PK part) did not change.

Record

Task 5

Reset all inputs. Set treatment start time to 2 and max drug efficacy to 1. Run the model.

Take a look at the remaining number of uninfected cells. Assume your goal was to prevent as many uninfected cells from getting infected as possible. For that purpose, would you want to try and double the half-life of the drug (i.e., reduce decay rate by a factor of 2) or reduce the concentration at which the drug has 50% action by half? Test both scenarios.

You should find that the for both changes, the number of uninfected cells is somewhat above 90K, and reducing the C50 value is slightly better (this of course depends on other model parameters, you can try to find a set of combinations of the other parameters at which the impact of these two hypothetical drug improvements switches).

Based on what you just found, which improved drug would you pick? If drug side-effects were an issue, would your decision change? Why or why not?

Record

Task 6

Let’s briefly consider a chronic situation. Reset all inputs, then set birth rate and death rate of uninfected cells to 10000 and 0.1. Run the simulation for 100 days.

You should get a steady state. Now, turn on drug with Emax = 1, treatment start at day 50, and a drug that’s given weekly at a dose of 1. Run the simulation.

You should see that once the drug is administered, the system is knocked out of its steady state, and settles back to another state. This is not a steady state, but instead a repeating pattern of virus increase and decrease. Make sure you understand why (take a look at the drug PK, you might need a log plot for that). Importantly, the drug can’t clear the infection.

Record

Task 7

Keep the settings as they were in the previous task. Play with the different PK parameters. Set C0 = 10. Not surprisingly, as the drug concentration increases, the virus load is reduced. Play around with the other PK parameters (e.g., dC and tinterval) to explore how changing those impacts the virus load.

Record

Task 8

Now, let’s revisit the PD parameters. Reset C0 = 1. Set C50 to 0.1. You should see a drop in the virus load compared to a value of 1. Next, let’s set C50 = 0.01. For this value, even small levels of drug are highly effective. You’ll see that for this scenario, the drug is strong enough to drive the virus towards zero/extinction, and the uninfected cells recover to almost their starting value.

Record

Task 9

Keep exploring. For both the acute and chronic infections, you can further explore how changes in PK (the model parameters that influence drug concentration) and PD (the model parameters that influence drug efficacy for a given concentration) lead to overall impact on the infection dynamics.

Record

Further Information

This app (and all others) are structured such that the Shiny part (the graphical interface you see and the server-side function that goes with it) calls an underlying R script (or several) which runs the simulation for the model of interest and returns the results.

For this app, the underlying function running the simulation is called simulate_pkpdmodel_ode. You can call them directly, without going through the shiny app. Use the help() command for more information on how to use the functions directly. If you go that route, you need to use the results returned from this function and produce useful output (such as a plot) yourself.

You can also download all simulator functions and modify them for your own purposes. Of course to modify these functions, you’ll need to do some coding.

For examples on using the simulators directly and how to modify them, read the package vignette by typing vignette('DSAIRM') into the R console.

If you want to learn a bit more about these kinds of models, see e.g. (Powers et al. 2003; Talal et al. 2006; Handel, Margolis, and Levin 2009; Canini et al. 2014).

References

Canini, Laetitia, Jessica M Conway, Alan S Perelson, and Fabrice Carrat. 2014. “Impact of Different Oseltamivir Regimens on Treating Influenza A Virus Infection and Resistance Emergence: Insights from a Modelling Study.” PLoS Computational Biology 10 (4): e1003568. https://doi.org/10.1371/journal.pcbi.1003568.
Handel, Andreas, Elisa Margolis, and Bruce R Levin. 2009. “Exploring the Role of the Immune Response in Preventing Antibiotic Resistance.” Journal of Theoretical Biology 256 (4): 655–62. https://doi.org/10.1016/j.jtbi.2008.10.025.
Powers, Kimberly A, Narendra M Dixit, Ruy M Ribeiro, Preeti Golia, Andrew H Talal, and Alan S Perelson. 2003. “Modeling Viral and Drug Kinetics: Hepatitis C Virus Treatment with Pegylated Interferon Alfa-2b.” Seminars in Liver Disease 23 Suppl 1: 13–18. https://doi.org/10.1055/s-2003-41630.
Talal, Andrew H, Ruy M Ribeiro, Kimberly A Powers, Michael Grace, Constance Cullen, Musaddeq Hussain, Marianthi Markatou, and Alan S Perelson. 2006. Pharmacodynamics of PEG-IFN Alpha Differentiate HIV/HCV Coinfected Sustained Virological Responders from Nonresponders. Hepatology 43 (5): 943–53.