Mathematica Modeling & Simulation Support for Clinical Trials

Experience and Examples from Phase Ib Studies of Hepatitis C Direct Antiviral Therapy

Robert B. (Nachbar^1)*, Jeffrey S. (Saltzman^1), Julie (Stone^2), Paul (Panorchan^2)

(( )^1)Applied Computer Science and Mathematics, Merck Research Laboratories, Rahway, NJ
(( )^2)Clinical PK/PD, Merck Research Laboratories, West Point, PA

*nachbar@merck.com

({{, , }})

Abstract

The design, execution, and analysis of a successful clinical trial requires many decisions to be made in a timely manner, often without all the hard data one would like to have to make informed choices. Modeling and simulation (M&S) is playing an expanding role to fill this gap. Mathematica offers many tools for M&S methods that can provide decisive information, visualization tools to effectively communicate the information, and documentation tools to adequately record, store, and retrieve the data and information. This presentation will describe many of the tasks undertaken with Mathematica to design a clinical trial, analyze the results, and modify the model in response to new information derived from the trial. We began by constructing a model for hepatitis C viral (HCV) infection based on the known biology of the infection process, and then converted it to a dynamical differential equation model. The model was combined with data from the literature for the antiviral interferon and used to predict the outcome of several Phase Ib trial designs for a Merck HCV NS3/4A protease inhibitor (PI). After the trial was completed and the data unblinded, a segmented regression model was constructed to allow nonlinear least squares fitting of viral load of all the patients simultaneously.  It  became obvious that the fitted model needed several changes to account for the pharmacokinetics of the PI and the emergence of resistance mutants that had not been anticipated prior to the trial. Interactive exploration of model and parameter values allowed us to quickly revise the kinetic model and refit it. Lastly, we used Monte Carlo methods to generate a large cohort of virtual patients, and visualization of the outcome of the simulated trial demonstrated that the model was able to predict the median and range of response to treatment in a dose-proportional way.

Introduction

Essentially, all models are wrong, but some are useful. — George E. P. Box [1]

The models should be based on the biology, chemistry, and/or physics of the problem

That is, we are not curve fitting

The model should be able to explain or predict a non-intuitive result

The models need only capture the essential features of the experiment (clinical trial) under investigation

Our models and the simulations derived therefrom are used to help make informed decisions

Is the treatment effective?

What is the best sampling strategy?

What is the most effective dose that is still safe?

This presentation will describe some of the work that was done to design and analyze a Phase 1b clinical trial for Merck's hepatitis C protease inhibitor MK-7009

Hepatitis C

World Health Organization [2]

Hepatitis C virus (HCV) is one of the viruses (A, B, C, D, and E), which together account for the vast majority of cases of viral hepatitis.

Globally, an estimated 170 million persons are chronically infected with HCV and 3 to 4 million persons are newly infected each year.

About 80% of newly infected patients progress to develop chronic infection.

Cirrhosis develops in about 10% to 20% of persons with chronic infection.

Liver cancer develops in 1% to 5% of persons with chronic infection over a period of 20 to 30 years.

HCV is spread primarily by direct contact with human blood. The major causes of HCV infection worldwide are use of unscreened blood transfusions, and re-use of needles and syringes that have not been adequately sterilized.

No vaccine is currently available to prevent hepatitis C.

Humans and chimpanzees are the only known species susceptible to infection, with both species developing similar disease.

An important feature of the virus is the relative mutability of its genome, which in turn is probably related to the high propensity (80%) of inducing chronic infection.

Hepatitis C Virus Infection

({{, (1)  virus binding and internalization
(2)  cytoplasmic release and uncoating
(3)  internal ribosomal entry site (IRES)-mediated translation and polyprotein processing
(4)  RNA replication
(5)  packaging and assembly
(6)  virion maturation and release.}})

A cartoon of the HCV structural and nonstructural proteins at the endoplasmic reticulum (ER) membrane. HCV RNA replication occurs in a specific, self-constructed membrane, the membranous web (MW). Positive and negative strand RNA are each generated within the MW respectively leading to viral production and sustaining the viral production machinery.

Model

Mechanism

We will consider a model after Neumann et al.[3] with 3 species: healthy hepatocytes (𝒯), virions (𝒱), and infected hepatocytes (ℐ).

\[KernelIcon];<<Chemistry`ChemicalKinetics`

\[KernelIcon];model={
𝒯+𝒱(Overscript[→, β])ℐ+𝒱  (* infection *),
ℐ(Overscript[→, (1-𝒻[t]ϵ)p])ℐ+𝒱  (* production, modulated by therapy *),
𝒱(Overscript[⇥, c])  (* virion clearance *),
ℐ(Overscript[⇥, δ])  (* infected cell death *),
𝒯(Overscript[→, ρ])2𝒯  (* regeneration *)
}

We can represent the model as a Petri net

\[KernelIcon];ReactionNetworkPlot[model,ImageSize->4 72,BaseStyle->{16}]

ReactionNetworkPlot uses GraphPlot

The graphical depiction of the model serves several purposes:

To visualize the interactions of the different parts of the model

To check the structure of the model

Model

Differential equations

\[KernelIcon];{vars,{odes,ics}}=KineticModel[model,t,InitialConditions->{\[LeftDoubleBracketingBar]𝒯\[RightDoubleBracketingBar][t0]==(Subscript[𝕋, 0]),\[LeftDoubleBracketingBar]ℐ\[RightDoubleBracketingBar][t0]==(Subscript[𝕀, 0]),\[LeftDoubleBracketingBar]𝒱\[RightDoubleBracketingBar][t0]==(Subscript[𝕍, 0])},InitialTime->t0]/.Amount->Identity;
Column[odes]

Since the volume of the mammalian liver is exquisitely controlled, we can impose the constraint

\[KernelIcon];totalHepatocyteSoln=Solve[𝒯[t]+ℐ[t]==(Subscript[𝔗, 0]),𝒯[t]]//First

Using this constraint on the total number of hepatocytes, we can eliminate either 𝒯[t] or ℐ[t].

\[KernelIcon];vars=DeleteCases[vars,𝒯];
odes=DeleteCases[odes,𝒯'[t]==_]/.totalHepatocyteSoln;
ics=DeleteCases[ics,𝒯[t0]==_];
Column[odes]

Model

Treatment naïve steady states

It will help at this point to introduce some nominal parameter values

\[KernelIcon];params={
(Subscript[𝔗, 0])->1.5×(10^7)(* cell/mL *),
p->6.6 (* (day^(-1)) *),
β->3.0×(10^(-7)) (* (((cell/mL))^(-1)) (day^(-1)) *),
δ->0.15 (* (day^(-1)) *),
c->6.2 (* (day^(-1)) *)};
paramNames=First/@params;

These values were established based on prior published work [3] and our own studies with chimpanzees.

The steady states of the system are found by solving for the variables with the (1^st) derivatives set to 0.

\[KernelIcon];SSeqn=odes/.(_)'[t]->0;
SSsoln=Solve[SSeqn,{ℐ[t],𝒱[t]}];

Let's first examine the situation without treatment:

\[KernelIcon];ssTable[SSsoln,params,{𝒻[t]->0,ϵ->0}]

We can determine the stability of these two steady states by examining the eigenvalues of the Jacobian of the system evaluated at each of the steady states:

\[KernelIcon];(J=Outer[D,Last/@odes/.𝒻[t]->0,{ℐ[t],𝒱[t]}])//MatrixForm

\[KernelIcon];λ=Eigenvalues[J]//Simplify

\[KernelIcon];λ/.SSsoln[[1]]/.𝒻[t]->0/.params//testStability

\[KernelIcon];λ/.SSsoln[[2]]/.𝒻[t]->0/.params//testStability

That is, in a healthy virtual individual an acute infection without treatment will always progress to a chronic infection, and an infected virtual individual will always remain infected.

This is in contrast to the real humans, where 25 % if acute infections spontaneously clear. Recall George Box's assertion about models!

Model

Steady states with treatment

Let's examine the situation with perfect treatment

\[KernelIcon];ssTable[SSsoln,params,{𝒻[t]->1,ϵ->1}]

\[KernelIcon];J=Outer[D,Last/@odes/.{𝒻[t]->1,ϵ->1},{ℐ[t],𝒱[t]}];
λ=Eigenvalues[J]//Simplify;

\[KernelIcon];λ/.SSsoln[[1]]/.params/.{𝒻[t]->1,ϵ->1}//testStability

\[KernelIcon];λ/.SSsoln[[2]]/.params/.{𝒻[t]->1,ϵ->1}//testStability

That is, with perfect treatment any infection will get cured, given enough time.

What if the treatment is less than perfect?

\[KernelIcon];ssTable[SSsoln,params,{𝒻[t]->1,ϵ->(97/100)}]

\[KernelIcon];J=Outer[D,Last/@odes/.{𝒻[t]->1,ϵ->(97/100)},{ℐ[t],𝒱[t]}];
λ=Eigenvalues[J]//Simplify;

\[KernelIcon];λ/.SSsoln[[1]]/.params/.{𝒻[t]->1,ϵ->(97/100)}//testStability

\[KernelIcon];λ/.SSsoln[[2]]/.params/.{𝒻[t]->1,ϵ->(97/100)}//testStability

So, a strong treatment (ϵ≥0.97) will lead to a cure.

What about a weak treatment?

\[KernelIcon];ssTable[SSsoln,params,{𝒻[t]->1,ϵ->(9/10)}]

An infected state is possible. Is it stable?

\[KernelIcon];J=Outer[D,Last/@odes/.{𝒻[t]->1,ϵ->(9/10)},{ℐ[t],𝒱[t]}];
λ=Eigenvalues[J]//Simplify;

\[KernelIcon];λ/.SSsoln[[1]]/.params/.{𝒻[t]->1,ϵ->(9/10)}//testStability

\[KernelIcon];λ/.SSsoln[[2]]/.params/.{𝒻[t]->1,ϵ->(9/10)}//testStability

That is, with a weak treatment, an infection will not be cured, regardless of the length of treatment.

Model

Cure boundary

We can find the boundary between success and failure by solving for ϵ when ℐ=0 and 𝒱=0.

\[KernelIcon];cbSoln=Solve[0=={ℐ[t],𝒱[t]}/.SSsoln[[2]]/.{𝒻[t]->1}//Thread,ϵ]

\[KernelIcon];cureBoundaryPlot[SSsoln[[2]]/.𝒻[t]->1,params,cbSoln]

Model

Viral dynamics under short-term therapy

We will assume that a patient enters a clinical trial with a steady state infection

\[KernelIcon];SSrules=SSsoln[[2]]/.{ℐ[t]->(Subscript[𝕀, 0]),𝒱[t]->(Subscript[𝕍, 0]),t->t0}

The current standard of care (SOC) for HCV is 48 weeks of pegylated Interferon α + Ribavirin with a 12 week follow-up period, so we can assume that on the time scale of weeks, the treatment is either on or off (𝒻[t]∈{0,1}). Here is one way to define a 7-day treatment

\[KernelIcon];treatment=𝒻[t_]:>(Piecewise[{{0, t≤0}, {1, 0<t≤7}, {0, 7<t}}]);

\[KernelIcon];vSoln=Table[NDSolve[{odes,ics/.SSrules/.t0->-7}/.treatment/.params,{ℐ[t],𝒱[t]},{t,-7,0,7,7+14}]//First,{ϵ,{0.9,0.95,0.99,0.999}}];
shortTermPlot[vSoln]

Trial Design

Parameter sensitivity

We want to design the clinical trial to provide the maximum amount of information

Therefore, we need to know how sensitive the model is wrt the parameters

We will use our package [4] for computing the derivatives of the variables wrt the parameters of the model.

\[KernelIcon];<<Chemistry`ODEDerivatives`

\[KernelIcon];?DefineODEDerivatives

\[KernelIcon];DefineODEDerivatives[vars, {odes, ics/.SSrules}/.treatment/.t0->-50, {t, -50, 25}, ℱ,Append[ paramNames,ϵ]]

\[KernelIcon];?ℱ

Here are some examples of the use of ℱ

\[KernelIcon];ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,ϵ][t]/.params

\[KernelIcon];%/.ϵ->0.999

\[KernelIcon];%/.t->-1

\[KernelIcon];{LogPlot[ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,ϵ][t]/.params/.ϵ->0.999,{t,-3,14},PlotRange->{Full,{(10^2),(10^8)}},ImageSize->4.5 72],Plot[D[ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,ϵ][t],c]/.params/.ϵ->0.999,{t,-3,14},PlotRange->{Full,{-3,1}(10^6)},ImageSize->4.5 72],Plot[c((D[ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,ϵ][t],c])/(ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,ϵ][t]))/.params/.ϵ->0.999,{t,-3,14},PlotRange->{Full,{-7,1}(10^0)},ImageSize->4.5 72]}//{{"𝒱","Sensitivity of 𝒱 wrt c","Normalized Sensitivity of 𝒱 wrt c"},#}&//Grid

Trial Design

Parameter sensitivity

\[KernelIcon];sensitivityPlot[]

Trial Design

Sampling times

We choose the sampling times where the viral load is most sensitive to the parameters.

Trial Design

5- vs. 7-Day study centers

Which days of the week are best for enrolling a patient?

({{, , }, {, , }})

For each enrollment day, the samples on the corresponding weekend days (red shading) were dropped from the synthetic data set (i.e., censored) and the model refit (blue curve). The 95 % confidence interval of the fit is shown as the dashed blue curves.

The only day of the week for enrollment that introduces large errors in the fit is Friday.

Preliminary Data Examination

MK-7009 Phase 1b clinical trial results [5]

({{, , }, {, , }})

Model Fitting

Parameter transformations and initial parameter estimates

Nonlinear optimization and regression generally need very good initial estimates. The values below came from Manipulate-ing the parameters by hand until good "fits" against the data above were obtained.

\[KernelIcon];initialParams=<<"initial parameter estimates.m"

We need a separate ϵ for each dose

\[KernelIcon];params=FilterRules[initialParams,Except[(Subscript[ϵ, _])]]

\[KernelIcon];epsilonParams=FilterRules[initialParams,(Subscript[ϵ, _])]

Because some of the (Subscript[ϵ, dose]) parameters are very close to 1, it will be advantageous to transform all the (Subscript[ϵ, dose]) to 1-(Subscript[ξ, dose]).

\[KernelIcon];xiParams=epsilonParams/.((Subscript[ϵ, d_])->x_):>((Subscript[ξ, d])->1-x)

\[KernelIcon];paramNames=First/@params

\[KernelIcon];xiParamNames=First/@xiParams

As (Subscript[ξ, 0]) is not adjustable, we define a rule for its evaluation

\[KernelIcon];(Subscript[ξ, 0.]):=1

Model Fitting

Regression function

We will use the special model function ℱ for constructing the regression model, but we redefine it to use ξ instead of ϵ.

\[KernelIcon];DefineODEDerivatives[vars, {odes, ics/.SSrules}/.ϵ->1-ξ/.treatment/.t0->-50, {t, -50, 25}, ℱ,Append[ paramNames,ξ]]
?ℱ

\[KernelIcon];LogPlot[ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,ξ][t]/.params/.ξ->(1-0.99)//Evaluate,{t,-7,21},PlotRange->{All,{1,(10^8)}}]

Model Fitting

Segmented model for multiple simultaneous responses

We're going to fit multiple simultaneous responses, one for each dose with its corresponding (Subscript[ξ, dose]). The other parameters will be common across all the responses.

We will construct the sum of squared differences of the experimental viral load values and those calculated by the model, and use FindMinimum to optimize the parameter values.

First load the clinical data for fitting

\[KernelIcon];fitData=<<"viral load data.m";
Dimensions[fitData]

\[KernelIcon];Join[fitData[[;;3]],{{"⋮","⋮","⋮"}},fitData[[-3;;]]]//TableForm[#,TableHeadings->{None,{"dose","time","HCV RNA"}}]&

Some observations were below the limit of quantitation (BLOQ), and are therefore censored.

\[KernelIcon];Cases[fitData,{__,_?(!NumberQ[#]&)}]//TableForm[#,TableHeadings->{None,{"dose","time","HCV RNA"}}]&

\[KernelIcon];{x,y}=With[{data=Cases[fitData,{__,_?NumberQ}]},{Thread[{dose,t}->#]&/@data[[All,{1,2}]],data[[All,3]]}];
First/@%

Because of the 6 orders of magnitude range in the response data, we will fit the logarithm of the viral load.

\[KernelIcon];minFunc=((((Log[y]-Log[(ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,(Subscript[ξ, dose])][t]/.x)]))^2)//Total)/(Length[y]-1);
minFunc/.(params//ShedUnits)/.xiParams

Special derivative rules are needed to coerce the ξ variants into the proper form.

\[KernelIcon];xiDerivativeRules={Derivative[der__][ℱ[v_]][T0_,p_,β_,δ_,c_,xi_][t_]:>(Derivative[der][ℱ[v]][T0,p,β,δ,c,ξ][t]/.ξ->xi)}

\[KernelIcon];grad=((∂_(#))minFunc&/@Join[paramNames,xiParamNames])/.xiDerivativeRules;
grad/.(params//ShedUnits)/.xiParams

Now, to keep the parameters positive during optimization, we transform them to logarithms, and make corresponding changes to the objective function and the gradient.

\[KernelIcon];?toLog

\[KernelIcon];{{logMinFunc,logGrad},logParams}=toLog[{minFunc,grad},Join[params,xiParams]];
logParamNames=First/@logParams

Reality check

\[KernelIcon];logMinFunc/.logParams

\[KernelIcon];logGrad/.logParams

Model Fitting

Least squares fitting of data

({optLogMinFunc,optLogParams}=Module[{step=0,t0,
purgeInterval=5},PurgeODEDerivatives[ℱ];t0=TimeUsed[];
FindMinimum[logMinFunc,List@@@logParams,Gradient->logGrad,Method->"Gradient",StepMonitor:>(step++;If[Mod[step,purgeInterval]==0,PurgeODEDerivatives[ℱ]])]])//Timing

\[KernelIcon];optLogParams=<<"optimized parameter estimates.m"

Now, back-transform the optimized parameters

\[KernelIcon];optParams=(optLogParams//fromLog)/.((Subscript[ξ, d_])->x_):>((Subscript[ϵ, d])->1-x)

\[KernelIcon];(Subscript[ϵ, 0.]):=0

Model Fitting

Fit vs. the data

({{, , }, {, , }})

Model Fitting

Assessment of Fit

Overall, the fit is good—captures the dose-proportional  viral load decline during treatment and return to baseline after treatment

However,

Observed initial viral load decline for low doses is greater than predicted by model

Observed viral load return to baseline after treatment stops is much more gradual than predicted by model

Observed viral load decline levels off after initial drop for all but the highest dose

Model Revision

Address initial decline ⟶ asymptotics

We can make the reasonable assumption that in the first day or so of the study, the fraction of the liver that is infected remains fairly constant, that is ℐ'=0 and ℐ=(Subscript[𝕀, 0]). We can then substitute (Subscript[𝕀, 0]) into 𝒱' which can then be solved exactly. This will also reduce the number of parameters to 2 [6,7].

\[KernelIcon];Vode=odes[[2]]/.𝒻[t]->1/.ℐ[t]->(Subscript[𝕀, 0])/.(SSrules/.𝒻[t0]->0)

\[KernelIcon];Vic=ics[[2]]/.(SSrules/.𝒻[t0]->0)/.t0->0

\[KernelIcon];soln=DSolve[{Vode,Vic},𝒱[t],t]//Simplify//First

This result can be further simplified by substituting the initial condition (Subscript[𝕍, 0])=((-c δ+p β (Subscript[𝔗, 0]))/(c β)), where (Subscript[𝕍, 0]) is the clinically determined initial viral load.

\[KernelIcon];soln=Solve[Join[soln/.Rule->Equal,{Vic/.𝒱[0]->(Subscript[𝕍, 0])}],𝒱[t]]//Simplify//First

Model Revision

Address initial decline ⟶ asymptotics

Example—fit initial viral load drop for 500 mg BID patients

\[KernelIcon];fitFunc=𝒱[t]/.soln

\[KernelIcon];fitParams={c->3.,ϵ->0.9}

Employ the usual transformations

\[KernelIcon];{logFitFunc,logParams}=toLog[fitFunc/.ϵ->1-ξ,fitParams/.(ϵ->x_):>(ξ->1-x)]

\[KernelIcon];𝒹ℴ𝓈ℯ=500.;
baselineData=Cases[fitData,{𝒹ℴ𝓈ℯ,t_,v_}/;t≤0:>v];
Length[%]

\[KernelIcon];V0=Mean[baselineData]//N

\[KernelIcon];fit=NonlinearModelFit[Cases[fitData,{𝒹ℴ𝓈ℯ,t_,v_}/;0<t≤2.25:>{t,Log[v]}],logFitFunc/.(Subscript[𝕍, 0])->V0//Log,List@@@logParams,t]

\[KernelIcon];optLogParams500=fit["BestFitParameters"]

\[KernelIcon];optParams500=fromLog[optLogParams500]/.(ξ->x_):>(ϵ->1-x)

\[KernelIcon];{fit["ANOVATable"],fit["ParameterTable"]}//List//Grid[#,RowAlignments->Top,ColumnSpacings->5]&

\[KernelIcon];asymptoticFitPlot[fitData,optParams500]

Model Revision

Address slow recovery ⟶ drug washout and dose-response efficacy

Another noticeable difference between the model and the data is the too rapid rebound of the viral load in the model after treatment is withdrawn. We had assumed that ϵ would immediately return to 0 at that time, but it is more likely that there is a slow washout of the inhibitor which induces a gradual decay of the efficacy. Therefore, a better way to handle the efficacy is through an exponential decay wash-out of the dose when treatment terminates coupled a dose-response model ((Subscript[E, max]) or Hill equation) for the efficacy (ϵ) [8].

\[KernelIcon];treatment=𝒻[t_]:> ((Piecewise[{{0, t<0}, {1, 0≤t<7}, {(E^(-(t-7) (Subscript[k, e]))), 7≤t}}]));doseResponse=ϵ 𝒻[t_]:>((((dose 𝒻[t]))^σ)/((((dose 𝒻[t]))^σ)+((Subscript[ED, 50])^σ)));

\[KernelIcon];params=Join[FilterRules[optParams,Except[(Subscript[ϵ, _])]],{(Subscript[k, e])->1,(Subscript[ED, 50])->20,σ->2}]

\[KernelIcon];washoutPlot[params]

And the effect on the viral dynamics

\[KernelIcon];vSoln=Table[PiecewiseNDSolve[{odes,ics/.SSrules/.t0->-7}/.doseResponse/.treatment/.params,{ℐ[t],𝒱[t]},{t,-7,7+14}]//First,{dose,{25.,75.,250.,500.,700.}}];
shortTermWashoutPlot[vSoln]

Model Revision

Address leveling off of decline ⟶ mutants

Another difference is that only at the 2 highest doses do the data show the steady decline during the second phase. 
At the lower doses, the viral load tends to bottom out and then rise. This is due to the presence of less fit mutant variants in the initial population of virions at about 0.01 %. These mutants are not as susceptible to the treatment as the wild-type virions, and as the wild-type population declines, the mutant population expands.

Here are 2 examples showing typical data for low and high doses. The blue curves are the total viral load, and the red curves are for the R155K mutant [9].

({{, }})

Model Revision

Model with wild-type and mutant species

We replace 𝒱 with 2 virion species (Subscript[𝒱, wt]) and (Subscript[𝒱, m]) for wild-type and mutant, respectively, and ℐ with (Subscript[ℐ, wt]) and (Subscript[ℐ, m]) [8]. We also add mutation steps to the model with the probability of mutation given as ζ. Finally, to handle the relative fitness of the mutant, we replace the production rate p with (Subscript[p, wt]) and (Subscript[p, m]).

\[KernelIcon];twoSpeciesModel={
𝒯+(Subscript[𝒱, wt])(Overscript[→, β])(Subscript[ℐ, wt])+(Subscript[𝒱, wt]),                                   𝒯+(Subscript[𝒱, m])(Overscript[→, β])(Subscript[ℐ, m])+(Subscript[𝒱, m]),                              (* infection *)
(Subscript[ℐ, wt])(Overscript[→, (1-𝒻[t](Subscript[ϵ, wt]))(1-ζ)(Subscript[p, wt])])(Subscript[ℐ, wt])+(Subscript[𝒱, wt]) ,     (Subscript[ℐ, m])(Overscript[→, (1-𝒻[t](Subscript[ϵ, m]))(1-ζ)(Subscript[p, m])])(Subscript[ℐ, m])+(Subscript[𝒱, m]) ,      (* production *)
(Subscript[𝒱, wt])(Overscript[⇥, c]) ,                                                          (Subscript[𝒱, m])(Overscript[⇥, c]) ,                                                   (* virion clearance *)
(Subscript[ℐ, wt])(Overscript[⇥, δ]) ,                                                           (Subscript[ℐ, m])(Overscript[⇥, δ]) ,                                                   (* infected cell death *)
(Subscript[ℐ, wt])(Overscript[→, (1-𝒻[t](Subscript[ϵ, wt]))ζ (Subscript[p, wt])])(Subscript[ℐ, wt])+(Subscript[𝒱, m]),               (Subscript[ℐ, m])(Overscript[→, (1-𝒻[t](Subscript[ϵ, m]))ζ (Subscript[p, m])])(Subscript[ℐ, m])+(Subscript[𝒱, wt])                (* mutation *),
𝒯(Overscript[→, ρ])2𝒯                                                                                                                             (* regeneration *)
}

Model Revision

Model with wild-type and mutant species

\[KernelIcon];ReactionNetworkPlot[twoSpeciesModel,ImageSize->8 72,BaseStyle->{14}]

Model Revision

Model with wild-type and mutant species

Differential equations

\[KernelIcon];{vars,{odes,ics}}=KineticModel[twoSpeciesModel,t,InitialConditions->{\[LeftDoubleBracketingBar]𝒯\[RightDoubleBracketingBar][t0]==(Subscript[𝕋, 0]),\[LeftDoubleBracketingBar](Subscript[ℐ, wt])\[RightDoubleBracketingBar][t0]==(Subscript[𝕀, wt,0]),\[LeftDoubleBracketingBar](Subscript[𝒱, wt])\[RightDoubleBracketingBar][t0]==(Subscript[𝕍, wt,0]),\[LeftDoubleBracketingBar](Subscript[ℐ, m])\[RightDoubleBracketingBar][t0]==(Subscript[𝕀, m,0]),\[LeftDoubleBracketingBar](Subscript[𝒱, m])\[RightDoubleBracketingBar][t0]==(Subscript[𝕍, m,0])},InitialTime->t0]/.Amount->Identity;
Column[odes]

\[KernelIcon];totalHepatocyteSoln=Solve[𝒯[t]+(Subscript[ℐ, wt])[t]+(Subscript[ℐ, m])[t]==(Subscript[𝔗, 0]),𝒯[t]]//First

\[KernelIcon];vars=DeleteCases[vars,𝒯];
odes=DeleteCases[odes,𝒯'[t]==_]/.totalHepatocyteSoln;
ics=DeleteCases[ics,𝒯[t0]==_];
Column[odes]

\[KernelIcon];doseResponse=(Subscript[ϵ, m_]) 𝒻[t_]:>((((dose 𝒻[t]))^σ)/((((dose 𝒻[t]))^σ)+((Subscript[ED, 50,m])^σ)));

Model Revision

Treatment naïve steady states

While we have done the following analysis symbolically, we will do it partially numerically here for clarity.

\[KernelIcon];params={
(Subscript[𝔗, 0])->2.3×(10^7)(* cell/mL *),
(Subscript[p, wt])->6.7 (* (day^(-1)) *),
(Subscript[p, m])->6.0 (* (day^(-1)) *),
β->6.7×(10^(-8)) (* (((cell/mL))^(-1)) (day^(-1)) *),
δ->0.97 (* (day^(-1)) *),
c->7.3 (* (day^(-1)) *),
ζ->0.0000125};

\[KernelIcon];SSeqn=odes/.(_)'[t]->0;
SSsoln=Solve[SSeqn/.params,{(Subscript[ℐ, m])[t],(Subscript[ℐ, wt])[t],(Subscript[𝒱, m])[t],(Subscript[𝒱, wt])[t]}]//Simplify;

\[KernelIcon];SSrules=With[{SSsoln=Solve[SSeqn/.{𝒻[t]->0,(Subscript[ϵ, wt])->0,(Subscript[ϵ, m])->0},{(Subscript[ℐ, m])[t],(Subscript[ℐ, wt])[t],(Subscript[𝒱, m])[t],(Subscript[𝒱, wt])[t]}]//Simplify},Select[SSsoln,MatchQ[#/.params,{(_->_?Positive)..}]&]]/.{(Subscript[ℐ, wt])[t]->(Subscript[𝕀, wt,0]),(Subscript[𝒱, wt])[t]->(Subscript[𝕍, wt,0]),(Subscript[ℐ, m])[t]->(Subscript[𝕀, m,0]),(Subscript[𝒱, m])[t]->(Subscript[𝕍, m,0]),t->t0}//First;
SSrules/.params

\[KernelIcon];ssTable[SSsoln,params,{𝒻[t]->0,(Subscript[ϵ, wt])->0,(Subscript[ϵ, m])->0}]

We can determine the stability of these steady states by examining the eigenvalues of the Jacobian of the system:

\[KernelIcon];J=Outer[D,Last/@odes/.𝒻[t]->0,{(Subscript[ℐ, m])[t],(Subscript[ℐ, wt])[t],(Subscript[𝒱, m])[t],(Subscript[𝒱, wt])[t]}]/.params;

\[KernelIcon];J/.SSsoln[[1]]/.𝒻[t]->0//Eigenvalues//testStability

\[KernelIcon];J/.SSsoln[[2]]/.𝒻[t]->0//Eigenvalues//testStability

\[KernelIcon];J/.SSsoln[[3]]/.𝒻[t]->0//Eigenvalues//testStability

As before, in a healthy virtual individual an acute infection without treatment will always progress to a chronic infection, and an infected individual will always remain infected.

Model Revision

Steady states with treatment

Let's examine the effects of different efficacies under treatment

\[KernelIcon];J=Outer[D,Last/@odes/.𝒻[t]->1,{(Subscript[ℐ, m])[t],(Subscript[ℐ, wt])[t],(Subscript[𝒱, m])[t],(Subscript[𝒱, wt])[t]}]/.params;

High efficacy against both species

\[KernelIcon];ssTable[SSsoln,params,{𝒻[t]->1,(Subscript[ϵ, wt])->0.9,(Subscript[ϵ, m])->0.9}]

\[KernelIcon];J/.SSsoln[[1]]/.{𝒻[t]->1,(Subscript[ϵ, wt])->0.9,(Subscript[ϵ, m])->0.9}//Eigenvalues//testStability

\[KernelIcon];J/.SSsoln[[2]]/.{𝒻[t]->1,(Subscript[ϵ, wt])->0.9,(Subscript[ϵ, m])->0.9}//Eigenvalues//testStability

\[KernelIcon];J/.SSsoln[[3]]/.{𝒻[t]->1,(Subscript[ϵ, wt])->0.9,(Subscript[ϵ, m])->0.9}//Eigenvalues//testStability

A cure can be effected

Low efficacy against both species

\[KernelIcon];ssTable[SSsoln,params,{𝒻[t]->1,(Subscript[ϵ, wt])->0.2,(Subscript[ϵ, m])->0.2}]

\[KernelIcon];J/.SSsoln[[1]]/.{𝒻[t]->1,(Subscript[ϵ, wt])->0.2,(Subscript[ϵ, m])->0.2}//Eigenvalues//testStability

\[KernelIcon];J/.SSsoln[[2]]/.{𝒻[t]->1,(Subscript[ϵ, wt])->0.2,(Subscript[ϵ, m])->0.2}//Eigenvalues//testStability

\[KernelIcon];J/.SSsoln[[3]]/.{𝒻[t]->1,(Subscript[ϵ, wt])->0.2,(Subscript[ϵ, m])->0.2}//Eigenvalues//testStability

A wild-type dominant infection will persist

High efficacy against wild-type, low efficacy against mutant

\[KernelIcon];ssTable[SSsoln,params,{𝒻[t]->1,(Subscript[ϵ, wt])->0.9,(Subscript[ϵ, m])->0.2}]

\[KernelIcon];J/.SSsoln[[1]]/.{𝒻[t]->1,(Subscript[ϵ, wt])->0.9,(Subscript[ϵ, m])->0.2}//Eigenvalues//testStability

\[KernelIcon];J/.SSsoln[[2]]/.{𝒻[t]->1,(Subscript[ϵ, wt])->0.9,(Subscript[ϵ, m])->0.2}//Eigenvalues//testStability

\[KernelIcon];J/.SSsoln[[3]]/.{𝒻[t]->1,(Subscript[ϵ, wt])->0.9,(Subscript[ϵ, m])->0.2}//Eigenvalues//testStability

A mutant dominant infection will persist

Model Revision

Viral load under treatment

\[KernelIcon];shortTermManipulate[]

[ this is a static image of the manipulate ]

Model Fitting

Fast forward

At this point we brought other modeling and simulation tools into play.

NONMEM [9] was used to fit a population-based model.

Fitting was done in stages:

First for the short-term asymptotic model to fit c and (Subscript[ED, 50,wt]) with (Subscript[ED, 50,m]) held at 500×(Subscript[ED, 50,wt]).

Then, holding c, (Subscript[ED, 50,wt]), and (Subscript[ED, 50,m]) fixed, δ and (Subscript[p, wt]) were fit with the full model.

Finally, (Subscript[ED, 50,m]) and (Subscript[p, m]) were adjusted by hand to fit the recovery of ((Subscript[𝒱, m])/((Subscript[𝒱, wt])+(Subscript[𝒱, m])))to baseline (literature data [10]) and the overall appearance of the visual predictive check [11] (next slide)

Model Results

Visual predictive check (Monte Carlo)

({{, , }, {, , }})

Heavy red line is median prediction from 500 virtual patients; heavy dashed red lines are 95 % prediction interval.

Lessons, conclusions, and follow-on work

Steady state stability analysis can be used to understand the model behavior under limiting cases of treatment.

Parameter sensitivity analysis can be used to locate optimal times at which to sample the viral load in patients.

Simulating censored observations showed that a less expensive trial design could be used without significant loss of information.

Parameter and data transformations may be needed to ensure numerical algorithms are well-behaved.

The initial model was wrong, but still useful, and it could be extended in a biologically relevant way.

Asymptotic analysis can be used to simplify the model (fewer parameters) and apply it to a limiting case of data to obtain better parameter estimates.

Monte Carlo simulation can be used to validate a model.

The 2-species model has since been used to fit data from a Phase IIa trial, and to help design another Phase Ib and a Phase IIb trial.

The techniques and algorithms developed for MK-7009 have been applied successfully to direct antiviral clinical trials.

Acknowledgments

Much the work described here could not have been accomplished without the support and input from others.

({{Erin Quirk, Jack Valentine, Nilou Mobashery}, {Jackie Gress, Anis Khan, Diana Brainard}, {David Olsen, Peng Sun, Steve Carroll}, {Steve Ludmerer, Nora Shire, Richard Barnard}, {Peggy Hwang, Hamish Wright, Matt Anderson}})

References

"Hepatitis C" World Health Organization, http://www.who.int/mediacentre/factsheets/fs164/en/.

Box, G. E. P. (http://en.wikiquote.org/wiki/George_E._P._Box).

Neumann, A. U., Lam, N. P., Dahari, H., Gretch, D. R., Wiley, T. E., Layden, T. J., and Perelson, A. "Hepatitis C Viral Dynamics in Vivo and the Antiviral Efficacy of Interferon-α Therapy", Science, 1998, 111, 222-333.

Nachbar, R. B., Pak, I., and Spann, A., "Kinetic Models in Biology and Chemistry" Wolfram Technology Conference 2006, Champaign, IL, 12-14 October 2006 (http://library.wolfram.com/infocenter/Conferences/6463/).

Lawitz, E., Sulkowski, M., Jacobson, I., Faruqui, S., Kraft, W.K., Maliakkal, B., Al-Ibrahim, M., Ghalib, R., Gordon, S., Kwo, P., Rockstroh, J., Miller, M.L., Hwang, P.M.T., Gress, J.O., and Quirk, E.K. "Safety, tolerability and antiviral activity of MK-7009, a novel inhibitor of the hepatitis C virus NS3/4a protease, in patients with chronic HCV genotype 1 infection" American Association For the Study Of Liver Diseases 59th Annual Meeting (the Liver Meeting), San Francisco, CA, 31 October-4 November 2008.

Panorchan, P., Nachbar, R.B., Saltzman, J.S., Khan, A.A., Olsen, D.B., Quirk, E.K., Barnard, R.J., Hwang, P.M.T., Wright, D.H., and Stone, J.A. "Evaluation of the dose-response relationship to short-term monotherapy with the HCV protease inhibitor, MK-7009" American Conference On Pharmacometrics, Mystic, CT, 4-7 October 2009.

Saltzman, J.S., Nachbar, R.B., Panorchan, P., Stone, J.A., and Khan, A.A. "Applying regular perturbation analysis to HCV viral load equations" Society For Industrial Applied Mathematics: Challenges And Frontiers Meeting, San Francisco, CA, 9-10 October 2009.

Stone, J.A., Panorchan, P., Nachbar, R.B., Khan, A.A., Olsen, D.B., Barnard, R.J., Hwang, P.M.T., and Saltzman, J.S. "Viral dynamics modeling of viral load and resistance data from short-term monotherapy with the HCV protease inhibitor, MK-7009" American Conference On Pharmacometrics, Mystic, CT, 4-7 October 2009.

Barnard, R., Feng, M., Ramos, G., Himmelberger, A., Adeniji-Adele, A., Greene, S., Miller, M., Graham, D., Ludmerer, S., Carroll, S., Olsen, D., Quirk, E., and Hazuda, D. "Analysis of anitviral resistance to MK-7009 in Protocol 004, a Phase 1b Monotherapy Study" (4^th) International Workshop on Hepatitis C: Resistance & New Compounds, Boston, MA 25-26 June 2009.

NONMEM, http://www.icondevsolutions.com/nonmem.htm.

Sarrazin, C., Kieffer, T. L., Bartels, D., Hanzelka, B., Müh, U., Welker, M., Wincheringer, D., Zhou, Y., Chu, H-M., Lin, C., Weegink, C., Reesink, H., Zeuzem, S., and Kwong, A. D. "Dynamic Hepatitis C Virus Genotypic and Phenotypic Changes in Patients Treated With the Protease Inhibitor Telaprevir" Gastroentrol. 2007, 132, 1767-1777.

Holford, N. http://www.fmhs.auckland.ac.nz/sms/pharmacology/holford/projects/vpc.aspx.

Initialization

Directory

NotebookDirectory[]//SetDirectory;

SetOptions

SetOptions[Plot,Axes->False,Frame->True,BaseStyle->{14,FontFamily->"Times"}];
SetOptions[LogPlot,Axes->False,Frame->True,BaseStyle->{14,FontFamily->"Times"}];
SetOptions[ListPlot,Axes->False,Frame->True,BaseStyle->{14,FontFamily->"Times"}];
SetOptions[ListLogPlot,Axes->False,Frame->True,BaseStyle->{14,FontFamily->"Times"}];

ssTable

ClearAll[ssTable]
ssTable[SSsoln_,params_,treatment_]/;MemberQ[treatment,ϵ->_]:=Module[{vars,ss,ssN,labels,dir,pn},vars=SSsoln[[1,All,1]];dir=(Piecewise[{{1, ϵ==1}, {-1, True}}])/.treatment;pn=First/@params;ss=Limit[vars/.SSsoln/.FilterRules[treatment,Except[{ϵ}]],Sequence@@FilterRules[treatment,{ϵ}],Direction->dir,Assumptions->Thread[pn>0]]//Chop;ssN=ss/.params;labels=((Piecewise[{{"Uninfected", MatchQ[#,{(0)..}]}, {"Infected", MatchQ[#,{__?Positive}]}, {"Unphysical", True}}]))&/@ssN;Grid[Prepend[Transpose[Prepend[ss//Simplify//Together,vars]],Prepend[labels,""]],ColumnAlignments->Left,ColumnSpacings->3,ColumnLines->{1,0},RowLines->{1,0}]]
ssTable[SSsoln_,params_,treatment_]/;MemberQ[treatment,(Subscript[ϵ, wt])->_]:=Module[{vars,ss,ssN,labels,dirWT,dirM,pn},vars=SSsoln[[1,All,1]]//Sort;dirWT=(Piecewise[{{1, (Subscript[ϵ, wt])==1}, {-1, True}}])/.treatment;dirM=(Piecewise[{{1, (Subscript[ϵ, m])==1}, {-1, True}}])/.treatment;pn=First/@params;ss=vars/.SSsoln//.treatment//Chop;ssN=ss/.params//Chop;labels=((Piecewise[{{"Uninfected", MatchQ[#,{(0)..}]}, {"Infected", MatchQ[#,{__?Positive}]}, {"Unphysical", True}}]))&/@ssN;Grid[Prepend[Transpose[Prepend[ss//Simplify//Together,vars]],Prepend[labels,""]],ColumnAlignments->Left,ColumnSpacings->3,ColumnLines->{1,0},RowLines->{1,0}]]

testStability

testStability[vals_]:=vals\[Implies]If[And@@Negative/@Chop[vals],Style["Stable",ColorData["HTML"]["RoyalBlue"],Bold],Style["Unstable",ColorData["HTML"]["Red"],Bold]]

cureBoundaryPlot

ClearAll[cureBoundaryPlot]
cureBoundaryPlot[SSsoln_,params_,cbSoln_]:=Module[{vars,cbSolnN=cbSoln/.params//First},vars=SSsoln[[All,1]];Plot[vars/.SSsoln/.params//Evaluate,{ϵ,0,1},Axes->False,Frame->True,FrameLabel->{"ϵ","concentration"},PlotLabel->"Cure Boundary",BaseStyle->{18,FontFamily->"Times"},GridLines->{{{ϵ,DotDashed}}/.cbSolnN,{{0,Dashed}}},PlotRange->{{0.79,1.01},{-2×(10^6),1.5 (10^7)}},ImageSize->6 72]]

shortTermPlot

shortTermPlot[vSoln_]:=Module[{loq=25.,zeroV=0.0000748502994011976`,t1=0,t2=7},LogPlot[𝒱[t]/.vSoln//Evaluate,{t,-3,21},Axes->False,Frame->True,FrameLabel->{"time (d)","HCV RNA (IU/mL)"},PlotLabel->"Viral Dynamics Under Short-Term Treatment",BaseStyle->{18,FontFamily->"Times"},PlotRange->{{-3,21},{(10^(-5)),(10^8)}},Prolog->{{GrayLevel[0.90],Rectangle[{t1,Log[(10^(-6))]},{t2,Log[(10^8)]}]},{Dashing[1],Line[{{-3,Log[zeroV]},{21,Log[zeroV]}}]},{Dashed,Line[{{-3,Log[loq]},{21,Log[loq]}}]}},Epilog->{Text[Style["LOQ",10],{-2.5,Log[loq]},{-1,-1.25}],Text[Style["cure",10],{-2.5,Log[zeroV]},{-1,-1.25}],Text[Style["ϵ = 0.9",15],{4,Log[(10^6)]},{-1,-1.25}],Text[Style["ϵ = 0.999",15],{1.5,Log[(10^4)]},{1,1.25}],{Text[Style["viral clearance",{14,Bold}],{2.75,Log[(10^7)]},{-1,-0.5}],{Thick,Arrow[{{2.5,Log[(10^7)]},{0.45,Log[(10^6.5)]}}]}},{Text[Style["hepatocyte clearance",14,Bold],{3,Log[(10^2.5)]},{0,0.5}],{Thick,Arrow[{{3,Log[(10^2.75)]},{4,Log[(10^3.75)]}}]}},{Inset[Style["reinstated production",14,Bold],{10.25,Log[(10^5.)]},{-1,0}],{Thick,Arrow[{{10,Log[(10^5.)]},{7.25,Log[(10^5.)]}}]}},{Text["reinfection",{10,Log[(10^6)]},{-1.1,0}],{Thick,Arrow[{{10,Log[(10^6)]},{7.5,Log[(10^6.75)]}}]}}},ImageSize->8 72]]

sensitivityPlot

sensitivityPlot[]:=Join[{With[{tmax=14,t0=-3,t2=7,loq=20,ymin=(10^0),ymax=(10^8)},LogPlot[ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,ϵ][t]/.params/.{{ϵ->0.9999`},{ϵ->0.999`},{ϵ->0.99`}}//Evaluate,{t,t0,tmax},Axes->False,Frame->True,FrameLabel->{None,StyleForm["𝒱 (IU (mL^(-1)))",FontSize->14]},PlotStyle->{{Hue[(16/24)]},{Hue[(8/24)]},{Hue[(0/24)]}},PlotRange->{{t0,tmax},{ymin,ymax}},MaxRecursion->15,PlotLabel->(StyleForm[#1,FontSize->14,FontColor->#2]&@@@Thread[{StringForm["ϵ = ``",ϵ]/.{{ϵ->0.9999`},{ϵ->0.999`},{ϵ->0.99`}},{Hue[(16/24)],Hue[(8/24)],Hue[(0/24)]}}]),Prolog->{GrayLevel[0.9],Rectangle[{0,Log[ymin]},{t2,Log[ymax]}]},Epilog->{{Text["LOQ",{t0,Log[loq]},{-1.5,-1},TextStyle->{FontSize->10}],Dashing[Medium],Line[{{t0,Log[loq]},{tmax,Log[loq]}}]}},ImageSize->4.5 72]]},MapThread[With[{tmax=14,t0=-3,t2=7,loq=20,ymin=#2,ymax=#3},Plot[#4(D[ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,ϵ][t],#]/ℱ[𝒱][(Subscript[𝔗, 0]),p,β,δ,c,ϵ][t])/.params/.{{ϵ->0.9999`},{ϵ->0.999`},{ϵ->0.99`}}//Evaluate,{t,t0,0,7,tmax},MaxRecursion->12,Axes->False,Frame->True,FrameLabel->{None,StyleForm[StringForm["((`1`)/𝒱)((∂𝒱)/(∂`1`))",#1],FontSize->14]},PlotLabel->StyleForm[#5,FontSize->14],PlotStyle->{{Hue[(16/24)]},{Hue[(8/24)]},{Hue[(0/24)]}},PlotRange->{{t0,tmax},{ymin,ymax}}(*,GridLines->{None,{0}}*),Prolog->{GrayLevel[0.9],Rectangle[{0,ymin},{t2,ymax}]},ImageSize->4.5 72]]&,{{p,β,δ,c,ϵ},{1.0,0.0,-1.3,-10,-1.1},{1.5,0.5,0.1,0.5,0.1},{p,β,δ,c,(ϵ/((1/(1-ϵ))1.5))},{"production (p)","infection (β)","hepatocyte clearance (δ)","viron clearance (c)","efficacy (ϵ)"}}]]//Partition[#,3,3,{1,1},{Graphics[{},AspectRatio->(1/GoldenRatio)]}]&//Grid[#,RowSpacings->2,ColumnSpacings->4]&//Grid[{{StyleForm["Sensitivity of the model to the parameters at three levels of efficacy",24,Bold,FontFamily->"Times"]},{#},{Graphics[{Text[StyleForm["time (d)",FontSize->14],{0,0}]},Frame->None,AspectRatio->(1/10)]}}]&

toLog, fromLog

ClearAll[toLog]
toLog::usage="toLog[model, params] returns a new model and parameter list with Log transformed parameters. params has the form {{(Subscript[par, i]), (Subscript[val, i])}..} or {(Subscript[par, i])->(Subscript[val, i])..}.";
toLog[model_,params:{{_,_}..}]:={model/.Thread[(First/@params)->(Exp/@log/@First/@params)],Thread[{(log/@First/@params),(Log/@Last/@params)}]}
toLog[model_,params:{__Rule}]:={model/.Thread[(First/@params)->(Exp/@log/@First/@params)],Thread[(log/@First/@params)->(Log/@Last/@params)]}
toLog[model_,allParams:{(𝒫[_,_]->_)..}]:=Module[{paramNames},paramNames=allParams[[All,1,2]]//Union;{model/.Thread[paramNames->Exp/@log/@paramNames],Thread[(𝒫[#1,log[#2]]&@@@First/@allParams)->(Log/@Last/@allParams)]}]

ClearAll[fromLog]
fromLog::usage="fromLog[logparams] backtransforms the Log transformed parameters. logparams has the form {{log[(Subscript[par, i])], (Subscript[val, i])}..} or {log[(Subscript[par, i])] -> (Subscript[val, i])}..}.";
fromLog[params:{{log[_],_}..}]:=Thread[{First/@params,(E^(Last/@params))}]/.log->Identity
fromLog[params:{(log[_]->_)..}]:=Thread[First/@params->(E^(Last/@params))]/.log->Identity
fromLog[params:{(𝒫[_,log[_]]->_)..}]:=Thread[First/@params->(E^(Last/@params))]/.log->Identity
fromLog[{}]:={}

asymptoticFitPlot

asymptoticFitPlot[data_,params_]:=Show[ListLogPlot[Cases[data,{𝒹ℴ𝓈ℯ,t_,v_}/;t≤2.25:>{t,v}],Axes->False,Frame->True,FrameLabel->{"time (d)","HCV RNA (IU/mL)"},PlotLabel->StringForm["`` mg BID",𝒹ℴ𝓈ℯ],PlotMarkers->Automatic,PlotRange->{{-1,3},{(10^1),(10^8)}},BaseStyle->18,ImageSize->6 72],LogPlot[V0,{t,-1,0},PlotStyle->{Thickness[Large],Red}],LogPlot[fitFunc/.(Subscript[𝕍, 0])->V0/.params//Evaluate,{t,0,3},PlotRange->Full,PlotStyle->{Thickness[Large],Red}]]

doseResponsePlot

doseResponsePlot[]:=Module[{ED50=15.,Emax=1.,σ={1,2,3}},LogLinearPlot[Emax((dose^σ)/((dose^σ)+(ED50^σ)))//Evaluate,{dose,(10^(-2)),(10^4)},Axes->False,Frame->True,FrameLabel->{"dose","response"},PlotLabel->"(Subscript[E, max]) Model",BaseStyle->{18,FontFamily->"Times"},PlotRange->All,GridLines->{{{ED50,DotDashed}},{{0.5Emax,Dashed}}},Epilog->{Inset[Style[#,24]&/@{"ϵ = (Subscript[E, max]) ((dose^σ)/((dose^σ) + ((Subscript[ED, 50])^σ)))","",StringForm["(Subscript[E, max]) = ``",Emax],StringForm["σ = ``",σ],StringForm["(Subscript[ED, 50]) = ``",ED50]}//Column,{Log[(10^(-0.5))],0.75Emax}]},ImageSize->8 72]]

washoutPlot

washoutPlot[params_]:=With[{ED50=15.,Emax=1.,σ={1,2,3},p={t1->0,t2->7},pd=(dose->{25.,75.,250.,500.,700.}//Thread//List//Thread),baseStyle={18,FontFamily->"Times"},imageSize=4 72},{Plot[dose 𝒻[t]/.treatment/.params/.p/.pd//Evaluate,{t,-3,0,15},Axes->False,Frame->True,FrameLabel->{"time (d)","dose (mg)"},PlotLabel->"Dose Washout",PlotRange->{Full,{-25,800}},ImageSize->imageSize,BaseStyle->baseStyle],LogLinearPlot[Emax((dose^σ)/((dose^σ)+(ED50^σ)))//Evaluate,{dose,(10^(-2)),(10^4)},Axes->False,Frame->True,FrameLabel->{"dose","response"},PlotLabel->"(Subscript[E, max]) Model",BaseStyle->baseStyle,PlotRange->All,GridLines->{{{ED50,DotDashed}},{{0.5Emax,Dashed}}},Epilog->{Inset[Style[#,12]&/@{"ϵ = (Subscript[E, max]) ((dose^σ)/((dose^σ) + ((Subscript[ED, 50])^σ)))",StringForm["(Subscript[E, max]) = ``",Emax],StringForm["σ = ``",σ],StringForm["(Subscript[ED, 50]) = ``",ED50]}//Column[#,Left,0.1]&,{Log[(10^(-2))],1.0Emax},{-1,1}]},ImageSize->imageSize],Plot[ϵ 𝒻[t]/.doseResponse/.treatment/.params/.p/.pd//Evaluate,{t,-3,0,15},Axes->False,Frame->True,FrameLabel->{"time (d)","ϵ"},PlotLabel->"Efficacy",PlotRange->{Full,{-0.05,1.05}},ImageSize->imageSize,BaseStyle->baseStyle]}]//List//Grid

shortTermWashoutPlot

shortTermWashoutPlot[vSoln_]:=Module[{loq=25.,zeroV=0.0000748502994011976`,t1=0,t2=7},Show[LogPlot[𝒱[t]/.vSoln//Evaluate,{t,-3,21},Axes->False,Frame->True,FrameLabel->{"time (d)","HCV RNA (IU/mL)"},PlotLabel->"Viral Dynamics Under Short-Term Treatment",BaseStyle->{18,FontFamily->"Times"},PlotRange->{{-3,21},{(10^(-6)),(10^8)}},Prolog->{{GrayLevel[0.90],Rectangle[{t1,Log[(10^(-6))]},{t2,Log[(10^8)]}]},{Dashing[1],Line[{{-3,Log[zeroV]},{21,Log[zeroV]}}]},{Dashed,Line[{{-3,Log[loq]},{21,Log[loq]}}]}},Epilog->{Text[Style["LOQ",10],{-2.5,Log[loq]},{-1,-1.25}],Text[Style["cure",10],{-2.5,Log[zeroV]},{-1,-1.25}],Text[Style["dose",10],{-2.5,Log[(10.^(-5.5))]},{-1,-1.25}],Text[Style["dose = 25",15],{3,Log[(10^5.5)]},{-1,-1.25}],Text[Style["dose = 700",15],{3,Log[(10.^2.5)]},{1,1.25}]}],Plot[Log[(10.^(-5.5))]+2ϵ 𝒻[t]/.doseResponse/.treatment/.{(Subscript[k, e])->2,(Subscript[ED, 50])->15,σ->2}/.dose->{25.,75.,250.,500.,700.}//Evaluate,{t,-3,21}],ImageSize->8 72]]

shortTermManipulate

shortTermManipulate[]:=Module[{soln,loq=25.,zeroV=0.0000748502994011976`,t1=0,t2=7,eqn,fixParams},fixParams=FilterRules[Join[params,{(Subscript[k, e])->1,(Subscript[ED, 50,wt])->15,(Subscript[ED, 50,m])->250,σ->1}],Except[{(Subscript[ED, 50,wt]),(Subscript[ED, 50,m]),(Subscript[k, e])}]];eqn={odes,ics/.SSrules/.t0->-7}/.doseResponse/.treatment/.fixParams;Manipulate[(soln=NDSolve[eqn/.{dose->p1,(Subscript[ED, 50,wt])->Exp[p2],(Subscript[ED, 50,m])->Exp[p3],(Subscript[k, e])->Exp[p4]},{(Subscript[ℐ, wt])[t],(Subscript[𝒱, wt])[t],(Subscript[ℐ, m])[t],(Subscript[𝒱, m])[t]},{t,-7,0,7,7+14}];
LogPlot[{(Subscript[𝒱, wt])[t],(Subscript[𝒱, m])[t],(Subscript[𝒱, wt])[t]+(Subscript[𝒱, m])[t]}/.soln//Evaluate,{t,-3,21},Axes->False,Frame->True,FrameLabel->{"time (d)","HCV RNA (IU/mL)"}(*,PlotLabel->"Viral Dynamics Under Short-Term Treatment"*),BaseStyle->{18,FontFamily->"Times"},PlotRange->{{-3,21},{(10^(-1)),(10^8)}},PlotStyle->{{Dashing[Medium],Thickness[Medium],Blue},{Dashing[Medium],Thickness[Medium],Red},{Thickness[Medium],Black}},Prolog->{{GrayLevel[0.90],Rectangle[{t1,Log[(10^(-6))]},{t2,Log[(10^8)]}]},{Dashing[1],Line[{{-3,Log[zeroV]},{21,Log[zeroV]}}]},{Dashed,Line[{{-3,Log[loq]},{21,Log[loq]}}]}},Epilog->{Text[Style["LOQ",10],{-2.5,Log[loq]},{-1,-1.25}]},ImageSize->8 72]),{{p1,250,"dose"},{25,75,250,500,700}},{{p2,10//Log,"(Subscript[ED, 50, wt])"},0.1//Log,50//Log},{{p3,250//Log,"(Subscript[ED, 50, m])"},10//Log,1000//Log},{{p4,1//Log,"(Subscript[k, e])"},0.5//Log,5//Log},SaveDefinitions->True]]

$Line

$Line=0;