LSD Pharmacokinetic Model

Kurt Bryan and SIMIODE

A notebook to start the analysis of the two-compartment pharmacokinetic LSD model in Section 6.5.1.

The Data: Sampling times, in hours

times = {1/12,1/4,1/2,1,2,4,8}

The concentration data (mg/liter) for each subject at the seven times above is shown below (with "0" for the missing data point at time t = 4 hours for subject 3).

concdata = {{11.1, 7.4, 6.3, 6.9, 5, 3.1, 0.8},{10.6, 7.6, 7, 4.8, 2.8, 2.5, 2},{8.7, 6.7, 5.9, 4.3, 4.4, 0, 0.3},{10.9, 8.2, 7.9, 6.6, 5.3, 3.8, 1.2},{6.4, 6.3, 5.1, 4.3, 3.4, 1.9, 0.7}};

A plot of the subject concentration data. First amalgamate data in to 5 lists of (times, concentration) data for each subject

concdata2 =Table[ Table[{times[[j]],concdata[[i,j]]},{j,1,7}],{i,1,5}];

ListPlot[concdata2,PlotMarkers->Automatic,PlotLegends->{"Subject 1","Subject 2", "Subject 3", "Subject 4", "Subject 5"}]

We use  cP(t) for the plasma concentration and cT(t) as the tissue concentration we have cP(0) = 12.27 for each subject and cT(0) = 0.

The ODE Model: The model of Section 6.5.1 (equation (6.66)) is

deP1 = 0.163*cP'[t] == 0.115*ka*cT[t] -0.163 *kb*cP[t] -0.163* ke*cP[t];
deT1 = 0.115*cT'[t] == 0.163*kb*cP[t] - 0.115*ka*cT[t];

where ka, kb, and ke are the constants to be determined from the concentration data.

As noted after equation (6.66), we can fill in the expressions

VP = 0.163*M;
VT = 0.115*M;

into the ODEs "deP1" and "deT1" above and find that "M" cancels. We are left with

deP =Assuming[M>0,Simplify[MultiplySides[deP1,1/M]]]
deT=Assuming[M>0,Simplify[MultiplySides[deT1,1/M]]]

These are equations (6.67) in the text.

We can have Mathematica solve these ODEs symbolically with initial data cP(0) = 12.27 and cT(0) = 0. The resulting solution is quite complicated, so the display is suppressed by putting a semicolon after DSolve (remove it if you want to see the solution.)

sol = DSolve[{deP,deT,cP[0]==12.27,cT[0]==0},{cP,cT},t];

Assign these to functions cPsol and cTsol:

cPsol=cP/. sol[[1]];
cTsol = cT/.sol[[1]];

The solution depends on ka, kb, and ke. We want to adjust these constants so that cP(t) fits the data.

Fitting the Data: Construct the least squares functional that depends on ka, kb, and ke. Omit missing data point for subject 3 at time 4 hours.

SS = Sum[Sum[(cPsol[times[[j]]]-concdata[[i,j]])^2,{i,1,5}],{j,1,7}];
SS = SS - (cPsol[times[[6]]]-concdata[[3,6]])^2;

Now we need to minimize SS as a function of ka, kb, and ke.  This can be done by setting partial derivatives equal to zero or by using Mathematica's "FindMinimium" command. It may help to provide an initial guess in which ka, kb, and ke are all positive.