The Art and Science of Model Translation

Robert B. Nachbar

Modeling & Simulation, Merck Research Laboratories

Abstract

Sharing information is one of the foundation stones of collaborative research. The exchange of mathematical and computational models between different software packages with perfect fidelity is required. Often times this task is difficult due to differences in the languages and the function libraries available to each, and at best it is tedious. Moreover, one does not necessarily need a full translation of an input file, as is done, for example, by the Unix command f2c for FORTRAN to C conversion. To meet these needs in our day-to-day modeling work, we have developed a set of practices and computational tools that facilitate the conversion of models from packages such as NONMEM, acslX, Berkeley Madonna, Jacobian, and MATLAB to Mathematica, and to a lesser extent from Mathematica to some of the others. Our methods make use of the symbolic computing capabilities of Mathematica, and demonstrate that code can, indeed, be treated as data. An additional benefit is that they permit one to rigorously prove that models independently implemented in different languages are identical, or not. The methods are based largely on string and expression replacement rules. We will show examples of the techniques applied to models for pharmacokinetics, disease progression, and ecology.

Motivation

Typical reasons for using someone else’s model:

The model is not in SBML or Modelica

You don't have the modeling application

Starting point for new research

Validate a model with new data

Explore model behavior with simulation and graphics

Extend the use beyond the capabilities of the native modeling application

Prove a model in X and your model in Mathematica are identical

The topic comes up frequently in on-line Mathematica forums

Programatically translating a model generally minimizes the potential for errors and produces a traceable record of the process

We will restrict the discussion to models composed of a system ordinary differential and algebraic equations with initial conditions

Strategic Questions & Tactical Issues

What are your needs?

Do you need to reproduce the whole computation, or only a part of it?

Is there potential for re-use of your translation code?

Does the other system/application employ a formal language?

Do you have the language specification?

How similar is it to Mathematica (Wolfram Language)?

We want to avoid writing a finite state automaton to parse the grammar

In most other systems the model description is procedural

Description of how to compute a step in the integration of the ODEs (NONMEM, MATLAB)

Order of operations can be important

Some modeling systems, such as acslX, handle a broader class of models, so our approach has been “just enough to get the job done right”

Fit for purpose, that is, don’t boil the ocean

Always with an eye toward reuse

Divide and conquer

Parameter assignments become a list of rules

Equation and initial condition statements become a list of equations

Integration information (e.g., method, initial and final time, etc.) are copied to NDSolve by hand

Caveat: The methods used here provide a very good starting point, but they are not comprehensive and have not been tested with all possible forms of input.

Methods for Translating to Mathematica

The “What”

lines of text in a file ⟹ expressions in a notebook

The “How”

Two levels of translation:

Translating individual statements, or lines of code

Parsing a whole input file

Two levels of manipulation:

Strings with StringReplace

Expressions with Replace, ReplaceAll, and/or ReplaceRepeated

Assumptions

The model input is syntactically and semantically correct

Translating Individual Statements

Input for most modeling systems is based on FORTRAN

Parentheses are used for grouping and for functions

Real numbers may be in floating point notation, e.g. 6.023E23

Variable and parameter names may or may not be case sensitive

Some characters for variable and parameter names are metacharacters in Mathematica and need to be replaced

Some variable and parameter names may be Mathematica function names and need to be changed

String patterns are used instead of regular expressions because they are similar to expression patterns

Strings are manipulated into a form that ToExpression can use

Use short-named helper function predicates, such as this one for parenthesis-free strings:

\[KernelIcon];?pfQ

Often the context of a string is important, as for "=" vs. "=="

\[KernelIcon];assignmentRule//InputForm

The default string pattern matching behavior takes the longest string, but sometimes the shortest string is needed, as in the argument list to a function

\[KernelIcon];First[mathRulesBM]//InputForm

\[KernelIcon];?bpQ

Some languages , such as C, ignore all whitespce, and some, such as Berkeley Madonna, do not; so handling whitespace can be tedious

\[KernelIcon];logicalTestRulesBM//First//InputForm

Expressions are further manipulated into the final form for the model

Parsing Whole Input Files

Different parts of the model may be collected into different sections of the model definition file (e.g., NONMEM, acslX)

Parameters may be initialized in one line and reassigned in a computation-dependent manner in following lines (e.g., NONMEM)

Sequential reassignment is often equivalent to use of a Piecewise function

Initial conditions may be assumed to be 0 unless specifically assigned (e.g., Berkeley Madonna)

Generally best to read the whole file with Import[filename, "Lines"] and work with the list of strings

Midazolam Pharmacokinetics—Berkeley Madonna

Berkeley Madonna [1] input files are a mix of binary and textual data, and it is generally best to save the model as a plain text file with the menu item File ▶ Save Equations As...

The input is a mathematical decription of the model, and the order of the statements is not important

The following example is from Sandra Visser’s doctoral thesis [4]

\[KernelIcon];(lines=Import["midazolam PKPD example Equations.txt","Lines"])//tfn

Preprocess the Strings

First, strip the comments and trim the lines

\[KernelIcon];lines=StringReplace[lines,";"~~___~~EndOfString->""]//StringTrim;

Next, let’s convert everything to lowercase so there will not be a collision with built-in Mathematica symbols

\[KernelIcon];lines=ToLowerCase[lines];

Are there any special characters we need to handle?

\[KernelIcon];StringCases[StringReplace[lines,NumberString->""],c_/;!StringMatchQ[c,LetterCharacter|WhitespaceCharacter]]//Apply[Union,#]&

Only the underscore (_) , which we will convert to a dollar sign ($)

\[KernelIcon];lines=StringReplace[#,"_"->"$"]&/@lines;

Extract the Useful Sections of the Input

The ODEs begin with "d/dt(var)="

\[KernelIcon];(odeLines=Cases[lines,line_/;StringMatchQ[line,"d/dt("~~(var__/;pfQ[var])~~")"~~WhitespaceCharacter...~~"="~~__]])//Column

The initial conditions begin with "init var="

\[KernelIcon];(iceLines=Cases[lines,line_/;StringMatchQ[line,"init"~~Whitespace~~__~~WhitespaceCharacter...~~"="~~__]])//Column

All the other lines with an equal sign can be considered as parameter assignments

\[KernelIcon];(paramLines=Cases[DeleteCases[lines,Alternatives@@Join[odeLines,iceLines]],line_/;StringMatchQ[line,__~~"="~~__]])//Column

Translate Strings to Valid Mathematica Input

ODEs and initial conditions

Replace the assignment ("=") with Equal ("=="). Note that we use the reserved names time and starttime

\[KernelIcon];(odeLines=StringReplace[odeLines,"d/dt("~~(var__/;pfQ[var])~~")"~~WhitespaceCharacter...~~"=":>var<>"'[time]=="])//Column

\[KernelIcon];(iceLines=StringReplace[iceLines,"init"~~Whitespace~~var__~~WhitespaceCharacter...~~"=":>var<>"[starttime]=="])//Column

The arithmetic in the ODEs and the initial conditions looks OK

Parameters

First convert the IF statement to a Piecewise function

\[KernelIcon];(paramLines=StringReplace[paramLines,ifThenElseRulesBM])//pwHL

Next, replace the assignment  ("=") with Rule  ("->")

\[KernelIcon];(paramLines=StringReplace[paramLines,assignmentRule])//Short

Convert Strings to Expressions

ODEs

\[KernelIcon];odes=ToExpression[odeLines];
vars=Cases[odes,(v_)'[time]:>v,∞]//Union

\[KernelIcon];varRules=Thread[vars->Through[vars[time]]]

\[KernelIcon];(odes=odes/.lhs_==rhs_:>lhs==(rhs/.varRules))//Column

Initial conditions

\[KernelIcon];(ices=ToExpression[iceLines])//Column

Parameters

\[KernelIcon];(params=ToExpression[paramLines]/.varRules)//Column

Prettying it up

At this point, one could make the variables and parameters more natural with rules such as these

{m1->(Subscript[M, 1]),m2->(Subscript[M, 2]),time->t,ec50->(Subscript[EC, 50]),midazolam$plasma$ngml->(Subscript[Midazolam, plasma]),…}

but it's not necessary:

Run the Model

\[KernelIcon];(eqns={odes,ices}//.params//Flatten)//Column

\[KernelIcon];intVar={time,starttime,stoptime}//.params

\[KernelIcon];soln=NDSolve[eqns,vars,intVar,MaxStepFraction->(1/50)]//First

Plot the Results

Model variables

Here we plot the ODE solution along with the infusion function that drives the system

\[KernelIcon];modelPlotBM[{midazolam$plasma$ngml,midazolam$tissue$ngml,Max[4infusion//Exp10,1]}//.params/.soln,{time,0,350}]

Here is the plot of the model variables from Berkeley Madonna

Biological effect

The biological effect is the “parameter” effect

\[KernelIcon];effectPlotBM[effect//.params/.soln,{time,0,350}]

Here is the plot of the biological effect from Berkeley Madonna

Details

//Show[#,ImageSize->8 72]&

//Show[#,ImageSize->8 72]&

HCV Viral Dynamics—NONMEM

NONMEM [2] input is comprised of sections, each with its own unique header string, and the sections can appear in any order

Constant parameters and initial conditions are in the $PK section

ODEs are in the $DES section

Adjustable parameters are in the $THETA section

The response variable is in the $ERROR section

Because the input is actually spliced into FORTRAN code at runtime and compiled, the input is actually a computational decription of the model, and the order of the statements is very important

Human intervention may be necessary to produce a faithful mathematical description of the model without excessive programming

We can, however, provide checks for many potential problems

NONMEM uses reserved names for the special variables in the model

ODE variables use the reserved names A(1), A(2), … in the ODEs and A_0(1), A_0(2), … in the initial conditions

Parameters to be optimized use the reserved names THETA(1), THETA(2), …; the initial estimates are in the section $THETA in ascending order

Random effects use the reserved names ETA(1), ETA(2), … and ERR(1), ERR(2), …. They can be set to 0 to model typical results

The response variable uses the reserved name Y

Two different forms can used for IF tests, one on a single line, the other over several lines

The following example is from Poland et al. [5]

\[KernelIcon];(lines=Import["HCV Viral Dynamics Model--Poland & Chittenden.ctl","Lines"])//tfn

Strip the comments, trim the lines, and delete blank lines

\[KernelIcon];lines=StringReplace[lines,";"~~___~~EndOfString->""]//StringTrim//DeleteCases[#,""]&;
Length[%]

Parse the Different Sections

Because whitespace is not used as a delimiter in the sections we need, it can be deleted, which will simplify the string patterns

The function getSectionNONMEM returns the lines from the input that are in the specified section

\[KernelIcon];(desLines=getSectionNONMEM[lines,"$DES"])//Column

Differential Equations

Recast the differential equation variables, their derivatives, and convert their assignmetns to equations

\[KernelIcon];desLines=StringReplace[desLines,specialVarRulesNONMEM];

Recast the remaining assignments to rules

\[KernelIcon];desLines=StringReplace[desLines,assignmentRule];

It will be easier to recast the multi-line long IF statement after the lines are coverted to expressions

Recast the single-line IF statements

\[KernelIcon];desLines=StringReplace[desLines,ifStatementRulesNONMEM];

Convert the strings to expressions

\[KernelIcon];des=ToExpression[StringReplace[desLines,mathRulesNONMEM]];

Details

With[{desLines=Style[#,Plain,Gray]&/@getSectionNONMEM[StringReplace[Import["HCV Viral Dynamics Model--Poland & Chittenden.ctl","Lines"],";"~~___~~EndOfString->""]//StringTrim//DeleteCases[#,""]&,"$DES"]},TabView[{"$DES lines"->(highlightLines[desLines,{1,2,3,4,5,6,7,8,9,10,11},Black]//Column),"precomp."->(highlightLines[desLines,{1,2,3,4,5,6,7,8},Blue]//Column),"ODEs"->(highlightLines[desLines,{-3,-2,-1},Green//Darker]//Column),"short IF"->(highlightLines[desLines,{2},Red]//Column),"long IF"->(highlightLines[desLines,{3,4,5,6,7},Orange]//Column)}]]

Conditional Reassignment

Often, a parameter is assigned a value contingent on some test, and a piecewise function is the appropriate mathematical construct, however FORTRAN does not have such a syntax. Frequently this is accomplished in NONMEM in the following manner

X = X0 - Y/Z
IF(X.LT.0) X=0
IF[X.GT.X0) X=X0

When these statements are converted to expressions the result is

x->({{\[Piecewise], ({{0, x<0}, {x0, x>x0}, {x0-(y/z), True}})}})

The difficulty here is that x is used in the test. Therefore, we have included a check for this condition and leave it up to the user to apply an appropriate correction

\[KernelIcon];{x->x0-y/z,If[x<0,x->0],If[x>x0,x->x0]}/.sequentialReassignmentRuleNONMEM

Recast a rule followed by If's containing a rule with the same lhs to a Piecewise function

\[KernelIcon];des=des//.sequentialReassignmentRuleNONMEM;

Recast muli-line IF statements into a Piecewise function

\[KernelIcon];des=des//.multilineIfRuleNONMEM;

Time Dependence of the Treatment

The model as specified is not defined for t<0, and the peg-IFN + RBV treatment is not time-dependent. So let’s fix these issues

\[KernelIcon];(des=des/.{t≤eot->0<t≤eot,pegr->({{\[Piecewise], ({{pegr, 0<t≤eotpr}, {0, True}})}})}/.{(effect:(effr|xpsr)->Piecewise[x_,def_]):>(effect->Piecewise[Prepend[x,{0,t≤0}],def])}//Quiet)//Column

ODEs in Final Form

Finally, substitute the rules into the equations

\[KernelIcon];(odes=Cases[des,_Equal]//.Cases[des,_Rule])//Column

Extract the differential equation variables from the ODEs

\[KernelIcon];vars=Cases[odes,(v_)'[t]==_:>v,∞]//Union

Parameters and Initial Conditions

Get the desired section

\[KernelIcon];(pkLines=getSectionNONMEM[lines,"$PK"])//Column

Recast the initial conditions, optimization parameters, and random effects as variables, and the initial condition assignments as equations

\[KernelIcon];pkLines=StringReplace[pkLines,specialVarRulesNONMEM];

Recaset the remaining assignments to rules

\[KernelIcon];pkLines=StringReplace[pkLines,assignmentRule];

Recast the IF tests

\[KernelIcon];pkLines=StringReplace[pkLines,ifStatementRulesNONMEM];

Convert the strings to expressions

\[KernelIcon];pk=ToExpression[StringReplace[pkLines,mathRulesNONMEM]];

Recast a rule followed by If’s containing a rule with the same lhs to a Piecewise function, and mulit-line IF statements to Piecewise functions

\[KernelIcon];pk=pk//.{sequentialReassignmentRuleNONMEM,multilineIfRuleNONMEM};

Set the random effects to 0

\[KernelIcon];pk=pk/.{eta[_]->0,err[_]->0};

Initial Conditions and Parameters in Final Form

Finally, extract the rules for the initial conditions, convert them to equations, and substitute the remaining rules into the equations

\[KernelIcon];(ices=Cases[pk,(a[i_][0]==rhs_)]//.(params=Cases[pk,_Rule]))//Column

Notice that the initial conditions for a[3] depends on a[2][0]

\[KernelIcon];(ices=With[{lhs=First/@ices,rhs=Last/@ices},Thread[lhs==(rhs/.(ices/.Equal->Rule))]])//Column

These initial conditions are for a treatment-free infected steady state, and as we have already defined the model to have no drug effect for t<0 we can safely modify the initial time to be less than 0

\[KernelIcon];(ices=ices/.a[i_][0]->a[i][-2])//Column

Optimization Parameters

\[KernelIcon];(thetaLines=getSectionNONMEM[lines,"$THETA"])//Column

\[KernelIcon];thetaLines=StringReplace[thetaLines,{"("->"{",")"->"}","fix"->""}];

\[KernelIcon];(optParams=Thread[Array[theta,Length[thetaLines]]->Replace[ToExpression[StringReplace[thetaLines,mathRulesNONMEM]],{{x_}:>x,{_,x_}:>x,{_,x_,_}:>x},{1}]])//Column

Experimental Parameters

STUDY, DOSE, and EOT come from the input data, so we need some typical values

\[KernelIcon];expParams={study->1,dose->200.,eot->6.,eotpr->0};

Run the Model

\[KernelIcon];(eqns={odes,ices}//.params/.optParams/.expParams//Flatten)//Column

\[KernelIcon];soln=NDSolve[eqns,vars,{t,-2,28},MaxStepFraction->(1/25)]//First;

Plot the response of the model

The model response Y for fitting is specified in the $ERROR section, which is treated much like the $PK section

\[KernelIcon];(errorLines=getSectionNONMEM[lines,"$ERROR"])//Column

Recast the ODE variables, optimization parameters, and random effects as variables, and the assignments as rules

\[KernelIcon];(errorLines=StringReplace[errorLines,Join[specialVarRulesNONMEM,assignmentRule]])//Column

Convert the strings to expressions and set the random effects to 0

\[KernelIcon];(error=ToExpression[StringReplace[errorLines,mathRulesNONMEM]]/.{eta[_]->0,err[_]->0})//Column

\[KernelIcon];plotModelNONMEM[{Exp10[y]}//.error//.params/.soln,{t,-2,28}]

Manipulate the Model

Let’s make this example a little more interesting

\[KernelIcon];

Module[{eqns={odes,ices}//.params/.optParamsM//Flatten,plotFuncs={Exp10[y]}//.error//.params,soln,tminW=-(2/7)},
Manipulate[soln=NDSolve[eqns/.{study->2,dose->d,eot->7eotW,eotpr->7eotprW},vars,{t,-2,7 36},MaxStepFraction->(1/50)]//First;(*Column[eqns[[;;3]]/.{study->2,dose->d,eot->7eotW,eotpr->7eotprW}]*)
LogPlot[plotFuncs/.soln/.t->7tW,{tW,tminW,28},Frame->True,FrameLabel->{"time (week)","HCV RNA (IU/mL)"},FrameTicks->{Automatic,Automatic(*,Automatic,Table[{(10^(4i)),i},{i,0,1,0.2}]*)},PlotStyle->{{Thickness[Large],Blue},{Thickness[Large],Green//Darker}},GridLines->{Automatic,(10^(Range[-5,6,2]))},PlotRange->{{tminW,tmaxW},{(10^(-5)),(10^7)}},Prolog->{{RGBColor[1,1,0.9,0.75],Rectangle[{0,(10^(-5))//Log},{eotprW,(10^7)//Log}]},{RGBColor[0.9,0.9,1.0,0.75],Rectangle[{0,(10^(-5))//Log},{eotW,(10^7)//Log}]},{Text["LOD",{tminW,10//Log},{-1.25,-1.25}],Dashing[Medium],Line[{{tminW,10//Log},{36,10//Log}}]},{Text["Cure",{tminW,(10^(-4))//Log},{-1.25,-1.25}],Dashing[None],Line[{{-1,(10^(-4))//Log},{36,(10^(-4))//Log}}]}},BaseStyle->{16},ImageSize->7.5 72],

{{d,200,"dose (mg)"},{0,50,100,200,400}},
{{eotW,1,"DAA Tx time (week)"},{1,4,12,24}},
{{eotprW,0,"peg-IFN + RBV Tx time (week)"},{0,12,24}},
Delimiter,
{{tmaxW,2,"plot range (week)"},{2,12,24,36}},

SaveDefinitions->True]]

Conclusions

Translating models to Mathematica is not overly difficult

One must pay particular attention to the details for computational models

The methods used here are extensible:

C—minimal coverage, needs

a push down automaton to handle nested {…}

use of embedded Module to handle scoping of variables

Jacobian

acslX—includes a push down automaton to parse nested sections in the file

Lastly, let’s recap the process with an example of another kind of model [6, 7]

\[KernelIcon];

Module[{size=2,show},show=Show[#,ImageSize->size 4 72]&;image1=picture1;image1b=EdgeDetect[image1,1.75,0.124]//ColorNegate;image1c=Module[{row1,row2,col1,col2,nrow,ncol},{ncol,nrow}=ImageDimensions[image1b];{row1,row2}={50,270};{col1,col2}={204,386};ImageTake[image1b,{row1,row2},{col1,col2}]//ImagePad[#,{{col1-1,ncol-col2},{nrow-row2,row1-1}},1]&];image1d=ReplacePixelValue[image1c,longListOfPoints->1];SlideView[Pane[#,size 4 72{1,Divide@@Reverse@ImageDimensions[image1]}]&/@show/@{
Graphics[{}],
image1,
Show[image1b,Graphics[{{Inset[Style["import model as lines",18 ,FontFamily->"Times"],{160,300},{1,0}]}}]],
Show[ImageCompose[image1c,{ImageSubtract[image1c,image1b]//ColorNegate//ColorReplace[#,Black->Red]&,0.75}],Graphics[{{Inset[Style["delete lines not needed",18 ,FontFamily->"Times"],{160,300},{1,0}],Arrow[{{92.45,290.2},{147.7,23.51}}],Arrow[{{92.45,290.2},{133.6,224}}]}}]],
image1c,
Show[ImageCompose[image1d,{ImageSubtract[image1d,image1c]//ColorNegate//ColorReplace[#,Black->Red]&,0.75}],Graphics[{{Inset[Style["manipulate lines and\ndo other fine tuning",18 ,FontFamily->"Times"],{160,300},{1,0}],Arrow[{{152.9,283.4},{241.9,210}}],Arrow[{{152.9,283.4},{287.8,134.4}}]}}]],
image1d,
Show[image1d,Graphics[{{Inset[Style["convert the lines to\nanother medium",18 ,FontFamily->"Times"],{160,300},{1,0}]}}]],
Sequence@@Module[{img,left=170,bottom=60,nrow1,ncol1,nrow2,ncol2},img=picture2;{ncol1,nrow1}=ImageDimensions[image1d];{ncol2,nrow2}=ImageDimensions[img];image2=ImagePad[img,{{left,ncol1-(left+ncol2)},{bottom,nrow1-(bottom+nrow2)}},1];Table[ImageCompose[image2,{image1d,f}],{f,1,0,-1/3}]//Rest],
Show[image2,Graphics[{{Inset[Style["experiment with\ndifferent colors",18 ,FontFamily->"Times"],{160,300},{1,0}],Arrow[{{165,300},{372,320}}],Arrow[{{165,300},{370,305}}]},{Inset[Style["rearrange some\nof the posture",18 ,FontFamily->"Times"],{160,125},{1,0}],Arrow[{{165,125},{302,271}}],Arrow[{{165,125},{353,154}}]}}]],
Sequence@@Module[{img,left=170,bottom=120,nrow1,ncol1,nrow2,ncol2},img=picture3;{ncol1,nrow1}=ImageDimensions[image1d];{ncol2,nrow2}=ImageDimensions[img];image3=ImagePad[img,{{left,ncol1-(left+ncol2)},{bottom,nrow1-(bottom+nrow2)}},1];MapThread[Show[#1,Graphics[{Inset[Style[#2,18 ,FontFamily->"Times"],{160,300},{1,0}]}]]&,{Table[ImageCompose[image3,{image2,f}],{f,1,0,-1/3}]//Rest,{"move to canvas","blend and apply the oils","Dazzle Your Audience\nWith the Model in\nMathematica!"}}]]
}]]

Acknowledgements

Sandra Visser (Merck) for the midazolam example

Matt Pillsbury (Merck) for helpful comments

References

http://www.berkeleymadonna.com

http://www.iconplc.com/technology/products/nonmem/

http://www.acslx.com

Visser, S. A. G., Huntjens, D. R. H., van der Graaf, P. H., Peletier, L. A., and Danhof, M. J. Pharm. Exper. Ther., 2003, 307, 765-775.

Poland, B., and Chittenden, J. http://www.certara.com/images/uploads/files/HCV_modeling_workshop.pdf.

Sovek, C. “The Art of Charles Benson” http://www.sovek.com/publications/articles/benson/.

Shutterstock.com, http://image.shutterstock.com/display_pic_with _logo/167056/96395081/stock-photo-wood-mannekin-waiting-for-meeting-sitting-on-chair-over-white-96395081.jpg.

Initialization

NotebookDirectory[]//SetDirectory

Needs["DifferentialEquations`PiecewiseNDSolve`"]

tfn:=TableForm[#,TableHeadings->Automatic]&

warnStyle[x_]:=Style[x,{Red,Bold}]

Exp10[x_]:=(10^x)

highlightLines[lines_,rows_,color_:Blue]:=Fold[ReplacePart[#,#2->(#1[[#2]]/.{Plain->Bold,Gray->color})]&,lines,rows]

optParamsM:=Join[FilterRules[optParams,Except[theta[6]]],{theta[6]->6.0}]

General

pfQ::usage="pfQ[string] returns True if string is free of parentheses.";
pfQ[str_String]:=StringFreeQ[str,"("|")"]
pfQ[_]:=False

bpQ::usage="bpQ[string] returns True if the parenteses in string are balanced.";
bpQ[str_String]:=Catch[Module[{n=0},Scan[Switch[#,"(",n++,")",n--;If[n<0,Message[bpQ::bktmop,str,"(","",""];Throw[False]]]&,Characters[str]];n==0]]
bpQ[_]:=False

assignmentRule={(a_/;a≠"=")~~"="~~(b_/;b≠"="):>a<>"->"<>b};

pwHL[{a___String,str_String/;StringMatchQ[str,__~~"Piecewise"~~__],b___String}]:={Skeleton[Length[{a}]],Style[str,Blue],Skeleton[Length[{b}]]}
pwHL[x_]:=x

Berkeley Madonna

logicalTestRulesBM={a:(")"|WhitespaceCharacter..)~~"and"~~b:("("|WhitespaceCharacter..):>a<>"&&"<>b,a:(")"|WhitespaceCharacter..)~~"or"~~b:("("|WhitespaceCharacter..):>a<>"||"<>b};

ifThenElseRulesBM={a:("="|WhitespaceCharacter..)~~"if"~~(b:("("|WhitespaceCharacter..)~~test__~~c:(")"|WhitespaceCharacter..))~~"then"~~(d:("("|WhitespaceCharacter..)~~trueResult__~~e:(")"|WhitespaceCharacter..))~~"else"~~(f:("("|WhitespaceCharacter..)~~falseResult__~~g:(")"|EndOfString)):>a<>"Piecewise[{{"<>d<>trueResult<>e<>", "<>b<>StringReplace[test,logicalTestRulesBM]<>c<>"}, {"<>f<>falseResult<>g<>", True}}]"};

cap[str_String]:=StringReplace[str,f_~~r___:>ToUpperCase[f]<>r]
cap[x_]:=x

mathFuncRule[f_String]:=With[{F=cap[f]},f~~"("~~ShortestMatch[arg__/;bpQ[arg]]~~")":>F<>"["<>arg<>"]"]
mathFuncRule[x_]:=x

mathRulesBM=mathFuncRule/@{"abs","log","log10","exp","sqrt"};

modelPlotBM[plotFuncs_,plotVar_]:=LogPlot[plotFuncs//Evaluate,plotVar,Frame->True,FrameLabel->{"time (min)","Midazolam conc. (ng/mL)",None,"Infusion (mg/mL)"},FrameTicks->{Automatic,Automatic,Automatic,Table[{(10^(4i)),i},{i,0,1,0.2}]},PlotStyle->{{Thickness[Large],Black},{Thickness[Large],Red},{Thickness[Large],Green//Darker}},GridLines->{Automatic,{1,10,100,(10^3),(10^4)}},PlotLegends->Placed[LineLegend[{"Plasma","Tissue","Infusion"},LegendFunction->Framed,Background->White],{{0.8,0.95},{0.5,1}}],PlotRange->{{0,350},{0.75,(10^4)}},Exclusions->None,PlotPoints->60,BaseStyle->{16},AspectRatio->(389/803),ImageSize->8.25 72]

effectPlotBM[plotFuncs_,plotVar_]:=Plot[plotFuncs//Evaluate,plotVar,Frame->True,FrameLabel->{"time (min)","Effect"},PlotStyle->{Thickness[Large],Black},GridLines->Automatic,PlotRange->{{0,350},{12,26}},Exclusions->None,PlotPoints->60,BaseStyle->{16},AspectRatio->(389/803),ImageSize->8 72]

NONMEM

getSectionNONMEM[lines:{__String},section_String]/;StringMatchQ[section,"$"~~__]:=Module[{pAll=First/@Position[lines,line_/;StringMatchQ[line,"$"~~__],{1},Heads->False],pSection=Position[lines,line_/;StringMatchQ[line,section~~___],Heads->False]},If[Last[pAll]≠Length[lines],pAll=Append[pAll,Length[lines]]];If[MatchQ[pSection,{}],Return[$Failed],pSection=pSection[[1,1]]];pSection=Cases[Partition[pAll,2,1],{pSection,_}]//First;Take[lines,pSection]//DeleteCases[#,line_/;StringMatchQ[line,"$"~~__]]&//ToLowerCase//StringReplace[#,Whitespace->""]&]

specialVarRulesNONMEM={
"dadt("~~i:DigitCharacter..~~")"~~"=":>"a["<>i<>"]'[t]==",
"a("~~i:DigitCharacter..~~")":>"a["<>i<>"][t]",
"a_0("~~i:DigitCharacter..~~")=":>"a["<>i<>"][0]==",
"a_0("~~i:DigitCharacter..~~")":>"a["<>i<>"][0]",
"theta("~~i:DigitCharacter..~~")":>"theta["<>i<>"]",
"eta("~~i:DigitCharacter..~~")":>"eta["<>i<>"]",
"err("~~i:DigitCharacter..~~")":>"err["<>i<>"]"
};

logicalTestRulesNONMEM={
".lt."->"<",
".le."->"<=",
".eq."->"==",
".ne."->"!=",
".ge."->">=",
".gt."->">",
".and."->"&&",
".or."->"||"
};

ifStatementRulesNONMEM={
"if("~~ShortestMatch[test__/;bpQ[test]]~~")then":>"If["<>StringReplace[test,logicalTestRulesNONMEM]<>", Null]",
"if("~~ShortestMatch[test__/;bpQ[test]]~~")"~~result__:>"If["<>StringReplace[test,logicalTestRulesNONMEM]<>", "<>result<>"]"
};

mathRulesNONMEM=Join[{
"**"->"^",
man:NumberString~~"e"~~s:("+"|"-")~~(exp:DigitCharacter..):>man<>"*10.^"<>s<>exp},
mathFuncRule/@{"abs","log","log10","exp","sqrt"}];

sequentialReassignmentRuleNONMEM={a___,lhs_->rhs_,Longest[ifs:(If[_,lhs_->_]..),\[Placeholder]],z___}:>{a,lhs->Module[{tests={ifs}[[All,1]],rhss={ifs}[[All,2,2]]},If[!FreeQ[tests,lhs],Print[warnStyle["Tests \!\(tests\) depend on \!\(lhs\)"]]];With[{body=Thread[{rhss,tests}]},Piecewise[body,rhs]]],z};

multilineIfRuleNONMEM={a___,If[test_,Null],lhs_->trueResult_,else,lhs_->falseResult_,endif,z___}:>{a,lhs->Piecewise[{{trueResult,test},{falseResult,True}}],z};

plotModelNONMEM[plotFuncs_,plotVar_]:=LogPlot[plotFuncs//Evaluate,{t,-2,28},Frame->True,FrameLabel->{"time (day)","HCV RNA (IU/mL)"},FrameTicks->{Automatic,Automatic(*,Automatic,Table[{(10^(4i)),i},{i,0,1,0.2}]*)},PlotStyle->{{Thickness[Large],Blue},{Thickness[Large],Green//Darker}},GridLines->Automatic,PlotRange->{{-2,28},{1,(10^7)}},Prolog->{RGBColor[0.9,0.9,1.0,0.75],Rectangle[{0,1//Log},{6,(10^7)//Log}]},BaseStyle->{16},ImageSize->7.5 72]

models

picture1=Import["stock-photo-wood-mannekin-waiting-for-meeting-sitting-on-chair-over-white-96395081.jpg"];
picture2=Import["Benson-p69-large.jpg"]//ImageResize[#,240]&;
picture3=Import["Benson-p68.jpg"]//ImageResize[#,240]&;
Thumbnail/@{picture1,picture2,picture3}//Row

picture1=Import["http://image.shutterstock.com/display_pic_with_logo/167056/96395081/stock-photo-wood-mannekin-waiting-for-meeting-sitting-on-chair-over-white-96395081.jpg"];
picture2=Import["http://www.sovek.com/images/articles/am89_10/Benson-p69-large.jpg"]//ImageResize[#,240]&;
picture3=Import["http://www.sovek.com/images/articles/am89_10/Benson-p68.jpg"]//ImageResize[#,240]&;
Thumbnail/@{picture1,picture2,picture3}//Row

longListOfPoints=#+{203,99}&/@{{3,86},{3,87},{3,88},{3,89},{3,90},{3,91},{3,92},{3,93},{3,94},{3,95},{3,96},{3,97},{3,98},{3,99},{4,83},{4,84},{4,85},{4,86},{4,99},{4,100},{4,101},{5,82},{5,83},{5,102},{6,81},{6,82},{6,103},{7,80},{7,81},{7,104},{8,80},{8,87},{8,88},{8,89},{8,95},{8,96},{8,97},{8,98},{8,104},{9,80},{9,85},{9,86},{9,89},{9,90},{9,94},{9,99},{9,100},{9,104},{9,105},{10,79},{10,84},{10,85},{10,90},{10,94},{10,100},{10,105},{11,79},{11,84},{11,90},{11,94},{11,100},{11,105},{12,79},{12,84},{12,90},{12,94},{12,101},{12,105},{13,79},{13,84},{13,90},{13,94},{13,101},{13,105},{14,79},{14,84},{14,90},{14,94},{14,100},{14,101},{14,105},{15,79},{15,84},{15,90},{15,94},{15,100},{15,105},{16,79},{16,84},{16,90},{16,94},{16,100},{16,105},{17,80},{17,85},{17,90},{17,94},{17,98},{17,99},{17,104},{18,80},{18,85},{18,90},{18,94},{18,95},{18,96},{18,97},{18,98},{18,104},{19,80},{19,81},{19,85},{19,90},{19,103},{19,104},{20,81},{20,82},{20,85},{20,90},{20,102},{20,103},{21,82},{21,83},{21,84},{21,85},{21,90},{21,101},{21,102},{22,83},{22,90},{22,99},{22,100},{22,101},{23,91},{23,92},{23,93},{23,94},{23,95},{23,96},{23,97},{23,98},{27,81},{27,82},{27,83},{27,84},{27,85},{27,86},{27,87},{27,88},{27,89},{27,90},{27,91},{27,92},{27,93},{27,94},{27,95},{27,96},{27,97},{27,98},{28,80},{28,99},{28,100},{28,101},{29,79},{29,80},{29,102},{29,103},{30,79},{30,103},{31,80},{31,104},{32,80},{32,81},{32,82},{32,83},{32,84},{32,85},{32,86},{32,87},{32,88},{32,89},{32,90},{32,91},{32,92},{32,93},{32,94},{32,95},{32,96},{32,97},{32,98},{32,104},{33,98},{33,99},{33,104},{33,105},{34,100},{34,105},{35,100},{35,105},{36,101},{36,105},{37,100},{37,105},{38,100},{38,105},{39,101},{39,102},{39,103},{39,104},{39,105},{41,82},{42,80},{42,81},{42,83},{42,84},{42,95},{42,96},{42,97},{42,98},{42,99},{42,100},{43,80},{43,84},{43,93},{43,94},{43,95},{43,101},{43,102},{44,80},{44,84},{44,93},{44,102},{44,103},{45,79},{45,84},{45,92},{45,103},{45,104},{46,79},{46,84},{46,91},{46,92},{46,104},{47,79},{47,84},{47,91},{47,96},{47,97},{47,98},{47,99},{47,104},{47,105},{48,79},{48,83},{48,91},{48,95},{48,96},{48,100},{48,105},{49,79},{49,83},{49,90},{49,95},{49,100},{49,105},{50,79},{50,83},{50,84},{50,90},{50,95},{50,101},{50,105},{51,79},{51,84},{51,90},{51,94},{51,101},{51,105},{52,79},{52,84},{52,90},{52,94},{52,101},{52,105},{53,79},{53,84},{53,90},{53,94},{53,101},{53,105},{54,80},{54,84},{54,85},{54,89},{54,94},{54,100},{54,105},{55,80},{55,85},{55,86},{55,87},{55,88},{55,89},{55,93},{55,94},{55,100},{55,105},{56,80},{56,81},{56,91},{56,93},{56,100},{56,105},{57,81},{57,93},{57,100},{57,104},{58,82},{58,92},{58,93},{58,101},{58,103},{58,104},{59,83},{59,84},{59,90},{59,91},{60,84},{60,85},{60,86},{60,87},{60,88},{66,110},{66,111},{67,112},{68,112},{68,113},{69,113},{70,113},{70,112},{71,112},{71,111},{71,104},{71,105},{72,105},{89,31},{89,32},{89,33},{89,34},{90,34},{91,34},{92,34},{93,34},{94,34},{94,33},{94,32},{94,31},{94,30},{94,29},{94,28},{93,28},{92,28},{92,29},{92,30},{93,31},{92,31},{91,31},{90,31}};