Project Help Notebook by Keith Stroyan

Help NoteBook Copyright 1998 by Academic Press, Inc. - All rights reserved.
Modify this Notebook to prepare your project report.

This notebook is to help you with the technical programming details of Project 35.  You must read the print project and develop the ideas from that.  This is only a shell document.

Project 35: Drug Dynamics & Pharmacokinetics

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

Initialization

Off[General::spell];
Off[General::spell1];
Off[InterpolatingFunction::dmval];
Off[ParametricPlot::pptr];
Off[Graphics::gprim];
Off[Graphics::gptn];
Off[NDSolve::ndsz];

Needs["Graphics`Colors`"]
	
colorQ[a_]:=(Head[a]===RGBColor);
grayQ[a_]:=(Head[a]===GrayLevel);
thickQ[a_]:=(Head[a]===Thickness);
absthickQ[a_]:=(Head[a]===AbsoluteThickness);
notruleQ[u_] := !(Head[u]===Rule);
ruleQ[u_] := (Head[u]===Rule);

Norm::usage = "Norm[expr] returns the Euclidean length of the
		vector, 'expr'.";
Norm[a_?VectorQ] := Sqrt[a.a];

tol = .5*10^-6;

Vector2D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics[Point[{0,0}]]];
	,						
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (.075 Norm[ntip-ntail]/Norm[ntip-ntail]) (ntip-ntail) ;
	up = 0.15 {-Punit[[2]],Punit[[1]]};
	arrow = Line[{ntail,ntip,ntip-Punit+up,ntip,ntip-Punit-up}];
	]];
Return[Graphics[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

		
uArrow[vec_,opts___]:= 
	Module[{a,b},
		If[Norm[N[vec]]<1,b = vec,b=vec/Norm[vec]];
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0}};
				Return[Vector2D[a,a+b,opts]]
				];

flow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
											
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

vectorFlow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
								
	arr[tt_] := Map[uArrow[{f[#[[1]],#[[2]]],g[#[[1]],#[[2]]]},Tail->#,Red]&,
	{x[t],y[t]}/.soln/.t->tt];
	
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0],arr[dt/100]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],arr[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

My Project Introduction

Write this last.  Begin with the basic question and state your main finding clearly, but with few technical details.  You might even modify the title of this section to make it more interesting.

Derivation of the Equations of Change

Variables & Parameters

Main Variables

Drug and Patient-specific Parameters

Amount vs. Concentration (HINT 35.3)

Concentration Drives Change, but Conservation of Matter Balances Things

Hints 35.4, 5, 6

The Final Result

Finally, we arrive at the result that once a drug has been introduced in a patient, or cB[0]  and  cT[0]  are given, they change according to the "laws:"

((Subscript[dc, B])/dt) =  -((Subscript[k, 1]) + (Subscript[k, 3])) (Subscript[c, B]) + (Subscript[k, 1]) (Subscript[c, T])

((Subscript[dc, T])/dt) =  (Subscript[k, 2]) (Subscript[c, B]) - (Subscript[k, 2]) (Subscript[c, T])

Our next section begins to explore the behavior of these equations.

Computer Experiments without Kidneys (HINT 35.7 ff.)

The Exact General Solution

The differential equations for blood and tissue concentration are linear equations in 2-D.  The text procedure of Chapter 24 could be applied to solvint them, but it would be pretty messy in the general case.  Mathematica doesn't mind the mess:

All Parameters Free, No Initial Conditions

Clear[f, g, cB, cT, t, k1, k2, k3, vB, vT];
f[cB_, cT_] := -(k1 + k3) cB + k1 cT;
g[cB_, cT_] := k2 cB  -  k2 cT;
DSolve[
      {cB'[t] == f[cB[t], cT[t]], 
        cT'[t] == g[cB[t], cT[t]]}, {cB[t], cT[t]}, t]

Mathematica can also solve the general case with a typical general initial contition as follows:

All Parameters Free, Initial Conditions

Clear[f, g, cB, cT, t, k1, k2, k3, vB, vT];
f[cB_, cT_] := -(k1 + k3) cB + k1 cT;
g[cB_, cT_] := k2 cB  -  k2 cT;
DSolve[
      {cB[0] == c0, cT[0] == 0, 
        cB'[t] == f[cB[t], cT[t]], 
        cT'[t] == g[cB[t], cT[t]]}, {cB[t], cT[t]}, t]

Mathematica can also solve specific cases exactly as follows:

Specific Parameters and Initial Conditions

Clear[f, g, cB, cT, t, k1, k2, k3, vB, vT];
f[cB_, cT_] := -(k1 + k3) cB + k1 cT;
g[cB_, cT_] := k2 cB  -  k2 cT;
k3 = 0;
k1 = 0.17;
vB = 2.10;
vT = 1.30;
k2  =  k1 vB/vT;
tfinal  =  10;
exactSoln = 
        DSolve[{cB[0] == 100, cT[0] == 0, 
            cB'[t] == f[cB[t], cT[t]], 
            cT'[t] == g[cB[t], cT[t]]}, {cB[t], cT[t]}, t];
cB[t_] = cB[t] /. exactSoln;
cT[t_] = cT[t] /. exactSoln;
Print["Blood concentration in red = ", cB[t]];
Plot[{cB[t], cT[t]}, {t, 0, tfinal}, PlotStyle -> {Red, Blue}];
      
Print["Tissue concentration in blue = ", cT[t]];

Mathematica can also solve specific cases numerically as follows:

Numerical Solution

Clear[f, g, cB, cT, t, k1, k2, k3, vB, vT];
f[cB_, cT_] := -(k1 + k3) cB + k1 cT;
g[cB_, cT_] := k2 cB  -  k2 cT;
k3 = 0;
k1 = 0.17;
vB = 2.10;
vT = 1.30;
k2  =  k1 vB/vT;
tfinal = 10;
approxSoln = 
        NDSolve[{cB[0] == 100, cT[0] == 0, 
            cB'[t] == f[cB[t], cT[t]], 
            cT'[t] == g[cB[t], cT[t]]}, {cB[t], cT[t]}, 
          {t, 0, tfinal}];
cB[t_] = cB[t] /. approxSoln;
cT[t_] = cT[t] /. approxSoln;
Print["Blood concentration in red."];
Plot[{cB[t], cT[t]}, {t, 0, tfinal}, PlotStyle -> {Red, Blue}];
      
Print["Tissue concentration in blue."];

HINT 35.9

Flow Solutions (HINT 35.10)

Clear[f, g, cB, cT, t, k1, k2, k3, vB, vT];
f[cB_, cT_] := -(k1 + k3) cB + k1 cT;
g[cB_, cT_] := k2 cB  -  k2 cT;
Print["Flow solution of the differential equations:"];
Print["dcB = ", f[cB, cT], " dt"];
Print["dcT = ", g[cB, cT], " dt"];
k3 = 0;
k1 = 0.17;
vB = 2.10;
vT = 1.30;
k2  =  k1 vB/vT;
initconds = Table[{cB, 0}, {cB, 10, 100, 10}];
tfinal = 10;
dt = 1;
Print["for the specific case:"];
Print["dcB = ", f[cB, cT], " dt"];
Print["dcT = ", g[cB, cT], " dt"];
soln = 
        flow2D[{f, g}, {cB, cT, t}, initconds, tfinal, dt, 
          AxesLabel -> {"cB", "cT"}, 
          PlotRange -> {{-1, 100}, {-1, 100}}];
Print["Explicit Plots of the various solutions."];
cB[t_] = cB[t] /. soln;
cT[t_] = cT[t] /. soln;
Plot[Evaluate[cB[t]], {t, 0, tfinal}, PlotStyle -> {Red}];
Plot[Evaluate[cT[t]], {t, 0, tfinal}, PlotStyle -> {Blue}];

HINT 35.12

HINT 35.13 Dahh, I actually can solve exactly.

Read the print project 35.13.  After you've seen it, you could have found it...

Summary of the patient with no kidneys

(Subscript[c, B])[t]  =  
      (Subscript[c, I]) ((Subscript[u, 1]) (e^(-h t)) + (Subscript[u, 2]))

(Subscript[c, T])[t]  =  
      (Subscript[c, I]) (-(Subscript[u, 2]) (e^(-h t)) + (Subscript[u, 2]))

where

(Subscript[u, 2]) = ((Subscript[v, B])/((Subscript[v, B]) + (Subscript[v, T])))	 and 	(Subscript[u, 1]) = ((Subscript[v, T])/((Subscript[v, B]) + (Subscript[v, T])))

COMPARE THESE FORMULAS WITH YOUR SPECIFIC COMPUTATIONS ABOVE

Flow Through Kidneys

Flow Solutions (HINT 35.14)

Clear[f, g, cB, cT, t, k1, k2, k3, vB, vT];
f[cB_, cT_] := -(k1 + k3) cB + k1 cT;
g[cB_, cT_] := k2 cB  -  k2 cT;
Print["Flow solution of the differential equations:"];
Print["dcB = ", f[cB, cT], " dt"];
Print["dcT = ", g[cB, cT], " dt"];
k3 = 0.03;
k1 = 0.17;
vB = 2.10;
vT = 1.30;
k2  =  k1 vB/vT;
initconds = Table[{cB, 0}, {cB, 10, 100, 10}];
tfinal = 10;
dt = 1;
Print["for the specific case:"];
Print["dcB = ", f[cB, cT], " dt"];
Print["dcT = ", g[cB, cT], " dt"];
flow2D[{f, g}, {cB, cT, t}, initconds, tfinal, dt, 
        AxesLabel -> {"cB", "cT"}, 
        PlotRange -> {{-1, 100}, {-1, 100}}];

Clear[f, g, cB, cT, t, k1, k2, k3, vB, vT];
f[cB_, cT_] := -(k1 + k3) cB + k1 cT;
g[cB_, cT_] := k2 cB  -  k2 cT;
k3 = 0.03;
k1 = 0.17;
vB = 2.10;
vT = 1.30;
k2  =  k1 vB/vT;
tfinal  =  20;
exactSoln = 
        DSolve[{cB[0] == 100, cT[0] == 0, 
            cB'[t] == f[cB[t], cT[t]], 
            cT'[t] == g[cB[t], cT[t]]}, {cB[t], cT[t]}, t];
cB[t_] = cB[t] /. exactSoln;
cT[t_] = cT[t] /. exactSoln;
Print["Blood concentration in red = ", cB[t]];
Plot[cB[t], {t, 0, tfinal}, PlotStyle -> {Red}];
Print["Tissue concentration in blue = ", cT[t]];
Plot[cT[t], {t, 0, tfinal}, PlotStyle -> {Blue}];

HINT 35.16 Equilibrium with Kidneys

Read and solve the print project hint.

A General Form for Exact Solutions HINT 35.17

Read the print project Hint 35.17 and solve it here.

Compare your textbook solution with the following:

All Parameters Free, No Initial Conditions

Clear[f, g, cB, cT, t, k1, k2, k3, vB, vT];
f[cB_, cT_] := -(k1 + k3) cB + k1 cT;
g[cB_, cT_] := k2 cB  -  k2 cT;
DSolve[
      {cB'[t] == f[cB[t], cT[t]], 
        cT'[t] == g[cB[t], cT[t]]}, {cB[t], cT[t]}, t]

Mathematica can also solve the general case with a typical general initial contition as follows:

More Generalities applied to Real Drugs and Patients

Hints 35.18, 19, 20, 21