SysModel applied for solving compartmental and physiological models

Author: Guillermo Sanchez:  guillerm@usal.es,   http://web.usal.es/~guillerm/

Update and additional material (for this tutorial is required Biokmod): http://web.usal.es/~guillerm/biokmod.htm

Here is introduced the compartmental analysis and it is shown how SysModel can be applied for solving compartmental model with constant coefficients.

We will need download these packages

Needs["Biokmod`SysModel`"]

Off[General::spell];
Off[General::spell1];
Needs["Graphics`Legend`"]
Needs["Graphics`Graphics`"]

Help about SysModel can be obtained as usual

?"Biokmod`SysModel`*"

Compartmental analysis has applications in clinical medicine, pharmacokinetics, internal dosimetry, nuclear medicine, ecosystem studies and chemical reaction kinetics.  It can be described as the analysis of a system in terms of compartments which separate the system into a finite number of component parts which are called compartments. Compartments interact through the exchange of species. Species may be a chemical substance, hormone, individuals in a population and so on.  A compartmental system is usually represented by a flow diagram or a block diagram.

We adopt the convention of  representing compartments with circles or rectangles. The flow into or out of the compartments is represented by arrows.  The (i^th) compartment  of a system of n compartments is labelled i and the size (amount or content) of the component in compartment  i as (Subscript[x, i])(t). The exchange between compartments, or between a compartment and the environment is labeled (Subscript[k, ij]), where i represents the flow from i to j. We can often write (Subscript[k, i])  if compartment i is not ambiguous with j. The environment is usually represented by "0" (zero), so (Subscript[k, i0]) is the fractional excretion coefficient from the i-th compartment to the outside environment.  The  input from the environment into the j-th compartment is called (Subscript[b, j])(t). Environment represents the processes that are outside the system. With regards to the environment, we only need to know the flow, (Subscript[b, j])(t),  into the system from the outside. The (Subscript[k, ij]) are called fractional transfer rate coefficients and they may be a function of different variables or constants.

The model of figure 1 consist in two compartments. There is an input (Subscript[b, i])(t)  to compartment i and an exchange of flow between compartment i and j, with transfer coefficients (Subscript[k, ij]) from compartment j to i  and  (Subscript[k, ji])  from compartment j to i. The model has an output of system from competent i with a transfer coefficient (Subscript[k, i0]).

The compartmental system is represented mathematically in terms of quantities such as the concentration or amount of a compound presents in the system. In the fig. 2  (Subscript[x, i])(t) and (Subscript[x, j])(t)  represent the quantities in compartment i and j at time t since the input (Subscript[b, i])(t) started. These quantities are called the state variables of the system and their evolution in the time are usually differential equations. So the system of the figure 2 can be described as it follow:

(Subscript[x, i])' (t) = -(Subscript[k, i0]) (Subscript[x, i])(t) - 
          (Subscript[k, ij]) (Subscript[x, i])(t) + (Subscript[k, ji]) (Subscript[x, j])(t) + b (t)
(Subscript[x, j])' (t) = 
        (Subscript[k, ij]) (Subscript[x, i])(t) - (Subscript[k, ji]) (Subscript[x, j])(t)

or in matrix notation

(({{(Subscript[x, i])' (t)}, {(Subscript[x, j])' (t)}})) =(({{-(Subscript[k, i0])-(Subscript[k, ij]) , (Subscript[k, ji])}, {(Subscript[k, ij]), -(Subscript[k, ji])}}))+(({{b(t)}, {0}}))

These models are System of Ordinary Differential Equation with Constant Coefficients (SODECC).

On many ocasions we are interested in models where the inputs (Subscript[b, i])(t) are acute in t  =  0, this is  (Subscript[b, i])(0) = (Subscript[b, i]) and (Subscript[b, i])(0) = 0 for t ≠ 0.  Another frequent situations are the ones  where the inputs (Subscript[b, i])(t) are continuous, this is  (Subscript[b, i])(t) = (Subscript[b, i]) for t  ≥ 0 or  (Subscript[b, i])(t) = (Subscript[b, i]) for  0 ≤ t  < T.

Figure 2 shows a more complex compartmental model

This model is represented by the following system of ODEs

\(TraditionalForm\`\(\(\(x\& . \_1\)(t)\)\(\ \)\(=\)\)\) (Subscript[b, 1])(
          t)  -  ((Subscript[k, 10])  +  (Subscript[k, 12])  +  (Subscript[k, 13])) (Subscript[x, 1])(
            t)   +  
        (Subscript[k, 21]) (Subscript[x, 2])(t)  +  (Subscript[k, 31]) (Subscript[x, 3])(
            t) 
\(TraditionalForm\`\(\(\(x\& . \_2\)(t)\)\(\ \)\(=\)\)\) (Subscript[b, 2])(t)  +  
        (Subscript[k, 12]) (Subscript[x, 1])(
            t)   -  ((Subscript[k, 20])  +  (Subscript[k, 21])  +  (Subscript[k, 23])) (Subscript[x, 2])(
            t)  +  (Subscript[k, 32]) (Subscript[x, 3])(t)
\(TraditionalForm\`\(\(\(x\& . \_3\)(t)\)\(\ \)\(=\)\)\) (Subscript[b, 3]) (t)  +  (Subscript[k, 13]) (Subscript[x, 1])(t)   +  
        (Subscript[k, 23]) (Subscript[x, 2])(
            t)  -  ((Subscript[k, 30])  +  (Subscript[k, 31])  +  (Subscript[k, 32])) (Subscript[x, 3])(
            t)

in matrix notation

(OverDot[x])(t) =A x + b(t), t ≥ 0

with

x\(TraditionalForm\`\((t)\) = \ {\(x\_1\)(t), \ \(x\_2\)(t), \ \(x\_3\)(t)}\^T\) 

b(t) = {(Subscript[b, 1])(t),  (Subscript[b, 2])(
                t), (Subscript[b, 3])(t)}^T

where  A is given by

A= (({{(Subscript[a, 11]), (Subscript[a, 12]), (Subscript[a, 13])}, {(Subscript[a, 21]), (Subscript[a, 22]), (Subscript[a, 23])}, {(Subscript[a, 31]), (Subscript[a, 32]), (Subscript[a, 33])}}))

with

(Subscript[a, 11]) = - ((Subscript[k, 10])  +  (Subscript[k, 12])  +  (Subscript[k, 13])) ; 
    (Subscript[a, 12])  =  (Subscript[k, 21]); (Subscript[a, 13]) =  (Subscript[k, 31]);
(Subscript[a, 21]) = (Subscript[k, 12]); 
    (Subscript[a, 22]) = - ((Subscript[k, 20])  +  (Subscript[k, 21])  +  (Subscript[k, 23]))  ;  
    (Subscript[a, 23]) = (Subscript[k, 32]);
(Subscript[a, 31]) = (Subscript[k, 13]);  (Subscript[a, 32]) =  (Subscript[k, 23]);  
    (Subscript[a, 33]) = - ((Subscript[k, 30])  +  (Subscript[k, 31])  +  (Subscript[k, 32]));

If we suppose that the substance introduced into the system is a radioactive isotope we mush considered the radioactive decay, it is given by a constant rate represented by λ. This constant is specific for each isotope. The decay constant can be interpreted by an equal flow going out of the system in each compartment. If we assume a radioactive decay in our example:

(Subscript[a, 11]) = - ((Subscript[k, 10])  +  (Subscript[k, 12])  +  (Subscript[k, 13]) +  λ) ; 
    (Subscript[a, 12])  =  (Subscript[k, 21]); (Subscript[a, 13]) =  (Subscript[k, 31]);
(Subscript[a, 21]) = (Subscript[k, 12]); 
    (Subscript[a, 22]) = - ((Subscript[k, 20])  +  (Subscript[k, 21])  +  (Subscript[k, 23]) +  λ)  ;  
    (Subscript[a, 23]) = (Subscript[k, 32]);
(Subscript[a, 31]) = (Subscript[k, 13]);  (Subscript[a, 32]) =  (Subscript[k, 23]);  
    (Subscript[a, 33]) = - ((Subscript[k, 30])  +  (Subscript[k, 31])  +  
            (Subscript[k, 32]) +  λ);

In some circumstances, a physiological model better a compartmental model is required, for example, when blood flow or other physiological parameter is measured. In those circumstances the transfer rate constants are associated with physiologically meaningful values that correspond to the measured physiological parameter or may be a function of them. However, if the physiological parameters are constants, a compartmental model and a physiological model are mathematically equivalent. But the form of the A matrix for the compartmental model compared with that for the physiological model will be different. For these cases we will directly give the values  (Subscript[a, ij]) to the A matrix instead of (Subscript[k, ij]).

The patterns that we have seen can be expanded to systems of n compartments or n state variables (in the case of physiological models). The equation for any compartment i is given by

(OverDot[x])(t) =A x + b(t), t ≥ 0
x(0) = (Subscript[x, 0])

where:

x\(TraditionalForm\`\((t)\) = \ {\(x\_1\)(t), \ \(x\_2\)(t), \  ... , \(x\_n\)(t)}\^T\) is a column vector and  (Subscript[x, i])(t) denote the amount or content of species in compartment i at time t.

A is a n×n is usually known as the compartmental matrix or system matrix.

b(t) = {(Subscript[b, 1])(t),  (Subscript[b, 2])(
                t),  ... , (Subscript[b, n])(t)}^T is a column vector where {(Subscript[b, i])(t)} is the input rate into compartment i from outside system.

x(0)  =  {(Subscript[x, 1])(0), (Subscript[x, 2])(0),  ..., (Subscript[x, n])(0)}\(\^T\)are the initial conditions, so (Subscript[x, i])(0), represents the amount or content of species in compartment i at time  t = 0.

SysModel.m provides a few functions to construct and solve Compartmental and Physiological models. Most of the examples are made using clearance rates, but also has been included some example of non linear models. The last example involved partial differential equations, it is solved using Laplace transform

The easiest example of compartmental model is the decay of a single radioactive isotope.

Let x be the number of isotope present in a sample at time t  and (Subscript[x, 0]) the number present at t = 0. The number of atoms decaying per unit of time is directly proportional present at time t, being  k the constant of proportionality. Then x(t) is given by the first order linear differential equation

d(x/d)t =  - 
          k x (t)

Eq.(2) can be evaluated using DSolve which gives the solution in analytical form. Initial conditions x[0] == x0 can have been added.

eq1 = DSolve[{x'[t]==-k x[t],x[0] == x0},x[t],t]

If  we know the values of k and (Subscript[x, 0]), for instant k  =  0.1 and (Subscript[(x ), 0]) =  100 with apropiated units, we can replace them in the previous expression to obtain x(t)

% gives the last result generated in the before output and expr /.  rules gives a rule of replacement  to be applied on expression expr.

% /. {k -> 0.1,  x0 ->  100}

In many occasions it is better build a Mathematica function where (Subscript[x, 0]) and k are parameters. They can be replaced with their values when ever it is necessary.

x[ t_,  x0_,  k_] = x[t] /. eq1[[1]]

Now we can replace k and x\(TraditionalForm\`\_0\) with their values or  k  =  0.1 and (Subscript[(x ), 0]) =  100

x[t, 0.1, 100]

Example 1.- Cobalt 60 is a radioactive isotope whose half-life is 5.3 year. It is used in medicine as a radiation source to kill cancer cells using gamma-emissions from cobalt 60. As it is radioactive its activity decreases with the time and the source must be changed periodically according to the expression above obtained x(t)  =  (Subscript[x, 0]) (e^(-k t)), where k is called the disintegration constant or transformation constant.

We are interested in evaluate the disintegration constant of \(TraditionalForm\`\^60\)Co. A half-life (T) of 5.3 year means that the radioactivity decrease to the half every 5.3 year, this is: if the activity is (Subscript[x, 0]) in t = 0 after  a time T the activity will be (Subscript[x, 0])/2.  Therefore we can calculate k doing x(t)  = (Subscript[x, 0])/2  and t =  T  in expression  x(t)  =  (Subscript[x, 0]) (e^(-k t)).

We use the function x[t, x0, k] before defined

eq =  Solve[ x[ T,  x0,  k]  ==  x0/2 ,  k]

Mathematica gives results in terms of inverse functions. In these cases it prints a warning message. We can suppress these warning messages

Off[Solve::"ifun"]; Off[InverseFunction::"ifun"];

This solution is applicable not only to the cobalt 60. It is a general relationship between k and T.

We replace k with its solution. We call it TransConst in order to be used later.

We use %%[[1]] to indicate to extract the result from two outputs before

TransConst  =  k /. %%[[1]]

As T is 5.3 year for Cobalt 60, and calling k = kCo60

kCo60  =  TransConst /. T -> 5.3

We can evaluate how the \(TraditionalForm\`\^60\)Co as function of time  decrease. So if we have 10 g of  \(TraditionalForm\`\^60\)Co in t  =  0, it will decrease as it follows.

We replace x0 with its value and using Plot we represent the evolution of Cobalt 60 in the time.

Plot[x[t, 10, kCo60],  {t, 0, 20}]

5 MBq of I-131 (Half-life 8.07 days) is needed for a diagnostic test.They are 3 days between shipment of radioiodine and its use in the test.How many MBq must be shipped?

The equation to resolve is x  =  (Subscript[x, 0]) (e^(-k t)) where x = 5 MBq and t = 3 days, k is de transformation constant of I-131 and (Subscript[x, 0]) is the variable to be calculated.

The transformation constant of I-131 is

kI131  = TransConst /. T -> 8.07

We replace t, k and x with their values.

Solve[ x[ 3 ,  x0, kI131 ]  ==   5 ,  
      x0]

Therefore should be shipped 6.47 MBq

In a experiment where a cell population is irradiated the survival of cells is given by eq.(1) with k = 1 (days^(-1)), and we are interested in calculating t for a specific value (Subscript[x, 1]) = 20% of initial population (Subscript[x, 0]).

We use Solve to obtain t.We can use the units as (days^(-1)). "%" indicating that % is a symbol and not the Math instruction %

Solve[x[t, 100 "%", 1 (days^(-1))] ==  20 "%", t] // 
      N

The following figure represents an isolated organ of constant volume with a liquid, with a dissolved drug,  flowing through the organ with a constant inflow and outflow rates v.

Let c(t) be concentration of the a drug in the organ. Let b(t) be the drug concentration in the liquid flowing into the organ. We assume a uniform and instantaneous mixing in the organ. The liquid enters with a v b(t)  rate and leave with a v x rate.  We can describe the process as follows: In the instant t the concentration is c(t), a later instant dt the variation of concentration in the organ will be: c+ dc   =  V c (drug in the organ in t) + c v dt (drug flowing organ) - q v dt (drug outflowing of the organ). Hence

V  c' (t) =  ν[
        b(t) - c(t)]

Let's consider the case for which the flowing concentration b(t) = b constant for and interval 0<t < T , and  b(t) =0 since  t > T,

eq2  =  V  c'[t] ==  ν (b - c[t]);

Let c(0) = (Subscript[c, 0]) be the concentration drug in the organ at t = 0.

DSolve[{eq2,  c[0] == c0},  c[t],  t]

We simplify c(t) and define  c1[t] in order to be used later

c1[t_] = c[t] /. FullSimplify[%][[1]]

If we know V, (Subscript[c, 0]) and b, we can replace them with their values in the above function, even plot the evolution of c(t). As example  we suppose that V = 5 L, (Subscript[c, 0]) = 0, ν = 0.2 L/s, b = 5 g/L

We plot c(t). We also use Timing that gives the time used in the computation.

Plot[
        c1[t] /. {V -> 5,  c0 -> 0,  ν ->  0.2,  
            b -> 5},  {t, 0, 100}] // Timing

Using Evaluate the time of computation is reduced, but cannot always we can apply this criterion (see apart. 1.9.1 of Mathematica Book).

Plot[
        Evaluate [
          c1[t] /. {V -> 5,  c0 -> 0,  ν ->  0.2,  
              b -> 5}],  {t, 0, 100}] // Timing

Suppose that the flow of the drug into the organ is cut in t  =  T.  Since that moment b(t)  =  0 but the liquid continue flowing with the same constant rates v. The process will be given by equation where t1  =  t - T

We do b  =  0 in eq2

V  c'[t1] ==  -ν c[t1];

The concentration in t1 = 0 (or t = T) will be with an initial concentration c[0] = c1

DSolve[{%,  c[0] == c1},  c[t1],  t1]

Replacing t1 = t-T and  c1 =  c1[T] is obtained the concentration c(t) in the organ  when t ≥ T

c2[t_, 
        T_] = c[t1] /. %[[1]] /. {t1 ->  t - T,  
          c1 ->  c1[T]}

Now we can build a joint function with c(t)  =  c1(t) when  0  ≤ t < T, and c(t)  =  c2(t) when t ≥ T

Remember that f[pattern/;condition]:= expr indicates that the expr is evaluated only when  condition is satisfied

c[t_, T_ /; t < T] := c1[t];
c[t_, T_ /; t ≥ T] :=  c2[t, T]

If T  =  50 the plot of c(t) is

Plot[
      c[t, 50] /. {V -> 5,  c0 -> 0,  ν ->  0.2,  
          b -> 5},  {t,  0, 100}]

We consider a bit more complicated situation: In eq.(2) we suppose that b(t)  =  A Exp(-k T) when  0  ≤  t  <  T, and b(t)  =  0 when t ≥ T .

The method is the same as we followed in the previous example. We have supposed that A= 1, k = 0.05

DSolve[{V  c'[t] ==  ν (A Exp[-k t] - c[t]),  
            c[0] == c0},  c[t],  t] // FullSimplify;

c3[t_] = c[t] /. %[[1]];

DSolve[{ V  c'[t1] ==  -ν c[t1],  
          c[0] == c1},  c[t1],  t1];

c4[t_, 
          T_] = c[t1] /. %[[1]] /. {t1 ->  t - T,  
            c1 ->  c3[T]};

c[t_, T_ /; t < T] := c3[t];
c[t_, T_ /; t ≥ T] :=  c4[t, T]

Plot[
      c[t, 30] /. {A ->  1,  k ->  0.05, V -> 5,  
          c0 -> 0,  ν ->  0.2,  b -> 5},  {t,  0, 
        100}]

Clear[x1, x2, eq,  x,  eq1, eq2, c, c1, c2, c3, c4,  
      kI131 ]

This Figure is represented an easy example of a two compartmental system (Anderson 1983). It shows the ingestion and metabolism of a drug. We have modeled with two compartments. It is supposed that the drug is taken orally flowing to the GI tract (Compartment 1) , then it is absorbed into the blood (Compartment 2) and finally eliminated.

We resolve the model of  fig. 1. Let (Subscript[x, 1])(t) and (Subscript[x, 2])(t), where t ≥ 0, is the mass of the drug in compartment 1 and 2, respectively.  Then we can describe the process as follow

d((Subscript[x, 1])/d)t =  
      b(t)  -  
        drug distribution rate from 1 to 2

d((Subscript[x, 2])/d)t =  
      inflow rate  -  
        outflow rate (elimination)

The previous equations are commonly called mass balance equation.  If it is assumed that rate of transference from each compartment i is proportional to the mass (or concentration) in this compartment, then (3) and (4) become

d((Subscript[x, 1])/d)t =  
      b(t)  - (Subscript[k, 12]) (Subscript[x, 1])

d((Subscript[x, 2])/d)t =  
      (Subscript[k, 12]) (Subscript[x, 1])  -  (Subscript[k, 20]) (Subscript[x, 1])

where (Subscript[k, 12]) and (Subscript[k, 20])  are the constants  (>0) of  proportionality  from 1 to 2 and from 2 to environment (elimination). This process is a simple case of first-order kinetics. 
Both ordinary differential equations (ODE) with appropriate initial conditions (Subscript[x, 1])(0) and (Subscript[(x ), 2])(0) constitute the compartmental metabolic model. In matrix-vector format the system of ordinary differential equation (SODE)  model is

(({{(Subscript[x, 1])' (t)}, {(Subscript[x, 2])' (t)}})) =(({{-(Subscript[k, 12]) , 0}, {(Subscript[k, 12]), (Subscript[k, 20])}}))+(({{b(t)}, {0}}))

This system can be solved with DSolve or using SystemDSolve. Here is assumed  k12= 0.3, and k20 = 0.05 and b(t) = {0, 0}, an as initial condition is taken {1,0}. The solution is plotted .

The first step is build the compartmental matrix

model1  =  
      CompartMatrix[2,  {{1, 2,  0.3},  {2, 0,  0.05}}]

{ x1[t_], x2[t_]}  =  {(Subscript[x, 1])[t], (Subscript[x, 2])[t]} /. 
        SystemDSolve[model1, {1, 0}, {0, 0}, t, t,  x]

Plot[{x1[t], x2[t]},  {t, 0, 100}, 
 
        PlotStyle -> {
{AbsoluteThickness[
                0.5]}, {AbsoluteThickness[1.5]}},  
        AxesLabel -> {"days",  "r"}, 
        PlotLegend -> {"x1", "x2"}, 
        LegendPosition -> {0.2,  -0.2},  
        LegendSize -> {0.3,  0.3}];

Clear[model1, x1, x2]

Below it is  shown a compartmental model with two inputs in compartment 1. One input represents a continuous flow  (Subscript[b, 1])(t) = (Subscript[b, 1]) for 0 < t < T and the other one consist in a single input  (bolus) (Subscript[b, 2])(t) =(Subscript[b, 0]) in t = 0 that it is equivalent to take as initial condition (Subscript[x, 1])(0) = (Subscript[b, 0]).

This process can be mathematically represented by

(Subscript[x, 1])' (t) = (Subscript[b, 1])(
          t)  -  ((Subscript[k, 10]) +  (Subscript[k, 12])) (Subscript[x, 1])(t)   +  
        (Subscript[k, 21]) (Subscript[x, 2])(t)  
(Subscript[x, 2])' (t)  =  
        (Subscript[k, 12]) (Subscript[x, 1])(t)   -   (Subscript[k, 21]) (Subscript[x, 2])(
              t) 
(Subscript[x, 1])(0) = (Subscript[b, 0]) 
(Subscript[b, 1])(t) = (Subscript[b, 1]) for 0 < t < T and  (Subscript[b, 1])(t) = 0 for t  ≥ T

Coefficients: (Subscript[k, 12])=2 (h^(-1)),  (Subscript[k, 21])=1 (h^(-1)), (Subscript[k, 10]) =0,5 (h^(-1)), 
Single input: (Subscript[b, 0]) = 100 mg
Constant input: (Subscript[b, 1])=10 mg/h until T = 10 h)

The first step is build the compartmental matrix

model2  = 
        CompartMatrix[
          2, {{1, 2, 2},  {2, 1, 1},  {1, 0, 0.5}}];

Here is the retention for  0 < t < T

ShowODE[model2,  {0, 0}, {10, 0},  t,  x1] // TableForm

sol1[t1_] = 
      SystemDSolve[model2,  {100,  0}, {10,  0},  t,  t1,  
        x1]

Here is the retention for  t  ≥ T.

\(ShowODE[
        model2, \ {\*""(Subscript[x1, 1])(T)"", \*""(Subscript[x1, 2])(T)""}, {0, 0}, \ t, \ \ x2] // TableForm\)

sol2[t_] = 
      AcuteInput[model2,  Map[Last, sol1[10]],  t,  x2]

We build a compound function that represents the retention divided in two intervals: one of them  when 0 < t < T and the other one t  ≥ T.

x[t_, T_] := 
        Map[Last, If[t < T,  sol1[t], sol2[t - T]]];

Here is the retention in compartment 1 and 2 f

Table[{t, {x[t, 10]}}, {t,  0, 20}] // Chop // 
      TableForm

Here is the retention in compartment 2

Plot[x[t, 10][[2]],  {t, 0, 12}]

Clear[model2,  sol1,  sol2, x1, x2,  x]

The figure shows a general two compartmental model.

Here is de compartmental matrix

model3 =CompartMatrix[
        2, {{1, 2,  k12},  {2, 1,  k21},  {1, 0,  k10},  {2, 0,  
            k20}}]

And here the SODE

ShowODE[model3,  {0, 0}, {b1, b2},  t,  x] // TableForm

In matrix notation can be represented as

where

x'(t) =(({{\(x\& . \_1\)}, {\(x\& . \_2\)}})) ; A =  (({{(Subscript[a, 11]), (Subscript[a, 12])}, {(Subscript[a, 21]), (Subscript[a, 22])}}));x =(({{(Subscript[x, 1])}, {(Subscript[x, 2])}})) ;  b(t) =(({{(Subscript[b, 1])}, {(Subscript[b, 2])}})) , t ≥ 0

with  (Subscript[a, 11])  = -((Subscript[k, 10]) + (Subscript[k, 12])) , (Subscript[a, 12 ]) =  (Subscript[k, 21]),  (Subscript[a, 21 ]) =  (Subscript[k, 12]),  (Subscript[a, 22]) = -((Subscript[k, 20]) + (Subscript[k, 21]))

In this case we can find a general  equation for any transfer coefficients  (with (Subscript[k, i]) ≠ (Subscript[k, j])).

\(\(\(\ \)\({x1[t1_], \ x2[t1_]} = \ {x\_1[t1], x\_2[t1]} /. 
          SystemDSolve[model3, \ {0, \ 0}, {b1, \ b2}, \ t, \ t1, \ x, \ 
            Method -> \ "DSolve"] (*\(//\)\(FullSimplify\)*) ;\)\)\)

The previous solution can be applied to a particular case substituting (Subscript[k, ij]) and the initial conditions (Subscript[X1, 0]) and (Subscript[X2, 0]) for their respective values, i.e.: Let's take (Subscript[k, 10 ]) =  0.05, (Subscript[k, 20]) = 0.3,  (Subscript[k, 12])  =  0.4, (Subscript[k, 21]) =  0.3,  (Subscript[X1, 0]) = 1, (Subscript[X2, 0]) = 1.

Plot[
      Evaluate[ {x1[t], x2[t]} /. {k10 -> 0.05,  k20 -> 0.03, 
            k12 -> 0.04, k21 -> 0.02,  b1 -> 1, 
            b2 -> 1 }],  {t,  0, 100}]

Clear[model3, x1, x2];

The analysis of model given above in the two variable cases can be inmediately generalizes to a system of with 3 or more compartments. In this section we show two models applied to represent drug administration process.  Other one is applied for people that intake lead in pollution enviromment. Finally it is presented the general tricompartmental model.

This figure show a tricompartmental model. In t  =  0 is assumed that one amount d is present in compartment 1 and 0 in the other compartments. A continuous liberation from a to 1 happen until that the initial amount d is completelly removed from a.

The kinetics of amount (Subscript[x, i])(t) will be given by the differential equations

\(x\& . \_a\)(t) = - (Subscript[k, a]) (Subscript[x, a])(t)
\(x\& . \_1\)(t) = (Subscript[k, a]) (Subscript[x, a])(t) - ((Subscript[k, 10])+ (Subscript[k, 12])) (Subscript[x, 1])(t)  + (Subscript[k, 21]) (Subscript[x, 2])(t)  
\(x\& . \_2\)(t) = (Subscript[k, 12]) (Subscript[x, 1])(t)  -  (Subscript[k, 21]) (Subscript[x, 2])(t)

We will add the initial conditions: (Subscript[x, a])(0) = d corresponds to an instantaneous "release"  d in t = 0. In the rest of the compartments we suppose that the concentration is null in t = 0, that is (Subscript[x, 1])(0) = 0; (Subscript[x, 2])(0) = 0;

model1 =  
      CompartMatrix[
        3, {{3, 1,  ka},  {1, 2, k12},  {2, 1,  k21},  {1, 0, 
            k10}}]

Thes equations are written as

eq1 = ShowODE[model1,  {d, 0, 0}, {0, 0, 0},  t,  x];

eq1 // TableForm // TraditionalForm

And here the solution is

{{x1[t_], x2[t_], xa[t_]}} = {(Subscript[x, 1])[t], (Subscript[x, 2])[t], (Subscript[x, 3])[t]} /. 
          DSolve[eq1, {(Subscript[x, 1])[t], (Subscript[x, 2])[t], (Subscript[x, 3])[t]}, t];

The solution is valid for any value of transfer coefficients. We apply for  the following values (Subscript[k, a])=2, (Subscript[k, 12]) = 2, (Subscript[k, 21])=1, (Subscript[k, 10]) = 0.5, d  =  100

Here is plot the concentration xa[t],  x1[t] and x2[t] for 0 ≤ t ≤ 50

Plot[
        Evaluate[{xa[t], x1[t], x2[t]} /. {d ->  100, ka -> 2, 
              k10 -> 0.5, k12 -> 2,  k21 -> 1}],  {t, 0, 
          30}, 
 
        PlotStyle -> {
{AbsoluteThickness[
                0.5]}, {AbsoluteThickness[1]}, {AbsoluteThickness[1.5],  
              AbsoluteDashing[{0.5, 1}]}},  
        AxesLabel -> {"days",  "r"}, 
        PlotLegend -> {"xa", "x1", "x2"}, 
        LegendPosition -> {0.2,  -0.2},  
        LegendSize -> {0.3,  0.3}];

Clear[xa, x1, x2, model1, eq1 ];

Other usual model in drug administration consist in a order-cero input to a constant rate and an absorption through a first order process. This process can be represented by a tricompartmental system as above section with we add and input  (Subscript[b, a])(t)  to the compartment as it show in the figure 8

We can represent the above model as

(Subscript[b, a])' (t)  =  - (Subscript[k, F]) 
(Subscript[x, a])' (t)  = 
        (Subscript[k, F])  -  (Subscript[k, a]) (Subscript[x, a])(t)
(Subscript[x, 1])' (t)  =  
        (Subscript[k, a]) (Subscript[x, a])(t)  -  ((Subscript[k, 10]) +  (Subscript[k, 12])) (Subscript[x, 1])(
              t)   +  (Subscript[k, 21]) (Subscript[x, 21])(t)  
(Subscript[x, 2])' (t)  =  
        (Subscript[k, 21]) (Subscript[x, 1])(t)   -   (Subscript[k, 12]) (Subscript[x, 2])(
              t)

Zero-order mean that the input is given by a constant rate (Subscript[b, a])(t)' = -(Subscript[k, F]), and therefore a constant input (Subscript[k, F]) to compartment a.   The total release amount is called d that is equivalent to taken as initial condition (Subscript[b, a])(0)  =  d.

(Subscript[b, a])[
        t_]  = (Subscript[b, a])[
          t]/.DSolve[{(Subscript[b, a])'[t] + (Subscript[k, F])  == 0, 
              (Subscript[b, a])[0] == d}, (Subscript[b, a])[t],  
            t][[1]]

We suppose that d = 100 and  (Subscript[k, F])=0.1  and plot (Subscript[b, a])[t]

Plot[
      Evaluate[(Subscript[b, a])[t] /. {d ->  100, (Subscript[k, F]) -> 0.1} ],  {t,  0, 
        1200}]

We will observe that from a certain value of t  (Subscript[b, a])(t) takes negative values.  This means the all d amount has been already transferred to a when  (Subscript[b, a])(t) = 0.  We solve this equation to find t that we will call T.

We used Solve to find T.

Solve[d - T (Subscript[k, F]) ==  0,  T]

Therefore we should be consider two period:

First period:  0<t≤T with T  =  d/(Subscript[k, F])

(Subscript[x, a])' (t)  =   (Subscript[k, F]) -  
          (Subscript[k, a]) (Subscript[x, a])(t)
(Subscript[x, 1])' (t)  =  
        (Subscript[k, a]) (Subscript[x, a])(t)  -  ((Subscript[k, 10]) +  (Subscript[k, 12])) (Subscript[x, 1])(
              t)   +  (Subscript[k, 21]) (Subscript[x, 2])(t)  
(Subscript[x, 2])' (t)  =  
        (Subscript[k, 12]) (Subscript[x, 1])(t)   -   (Subscript[k, 21]) (Subscript[x, 2])(
              t)

Initial condition: (Subscript[x, a])(0) = (Subscript[x, 1])(0) =(Subscript[x, 2])(0) = 0

Second period: t > T. We call t' = t - T

(Subscript[x, a])' (t')  = - (Subscript[k, a]) (Subscript[x, a])(
            t')
(Subscript[x, 1])' (t')  =  
        (Subscript[k, a]) (Subscript[x, a])(t')  -  ((Subscript[k, 10]) +  (Subscript[k, 12])) (Subscript[x, 1])(
              t')   +  (Subscript[k, 21]) (Subscript[x, 2])(t')  
(Subscript[x, 2])' (t')  =  
        (Subscript[k, 12]) (Subscript[x, 1])(t')   -   (Subscript[k, 21]) (Subscript[x, 2])(
              t')

Initial condition: (Subscript[x, a])(0) = (Subscript[x, a])(T);  (Subscript[x, 1])(0) =(Subscript[x, 1])(T); (Subscript[x, 2])(0) = (Subscript[x, 2])(T), where  (Subscript[x, a])(T), (Subscript[x, 1])(T), y (Subscript[x, 2])(T) is given by eq. 13 when t = T.

Now, we give the values of transfer coefficients. In our example they are  (Subscript[k, F])=0.1, (Subscript[k, a])= 2, (Subscript[k, 10])=0.5, (Subscript[k, 12])= 2, (Subscript[k, 21])= 1

We use k with /: (TagSet) to  make easier clear all coefficient values later. We will can delete its writing Clear[k] in place of write every (Subscript[k, i]) . Also we write (Subscript[k, F]) = 1/10 and (Subscript[k, 10]) = 1/2 using rational coefficients better than (Subscript[k, F]) = 0.1 and (Subscript[k, 10]) = 0.5 because with rational coefficients the computation is made using infinite precision.

k /: (Subscript[k, F]) = 1/10 ;k /: (Subscript[k, a]) = 2 ;k /: 
        (Subscript[k, 10]) = 1/2  ;k /: (Subscript[k, 21]) = 1 ; k /: 
        (Subscript[k, 12]) = 2 ;

Here is the solution to eq13

eq1  =  {(Subscript[x, a])'[t] ==   (Subscript[k, F]) -  (Subscript[k, a]) (Subscript[x, a])[t], x1'[
              t] == 
            (Subscript[k, a]) (Subscript[x, a])[t] -  ((Subscript[k, 10]) + (Subscript[k, 12])) x1[t] + 
              (Subscript[k, 21]) x2[t], x2'[t] == 
            (Subscript[k, 12]) x1[t] - (Subscript[k, 21]) x2[t]};

{x1a[t_], x11[t_], 
        x12[t_]} = {(Subscript[x, a])[t], x1[t], 
            x2[t]} /. DSolve[
              Join[eq1, {(Subscript[x, a])[0] == 0, x1[0] == 0, 
                  x2[0] == 0}], {(Subscript[x, a])[t], x1[t], x2[t]} ,  
              t][[1]] // FullSimplify

Plot[{x1a[t], x11[t], x12[t]},  {t, 0, 50}, 
      PlotRange -> All]

Here is the solution to eq14. We have to bear in mind that xa[0] =  (Subscript[x, a]), x1[0] =  (Subscript[x, 1]) and x2[0] =  (Subscript[x, 2]) correspond to the retention in the compartment a, 1 and 2  estimated with eq13 for t  =  T, that is  (Subscript[x, a]) =  x1a[T],  (Subscript[x, 1]) =  x11[T]and (Subscript[x, 2]) =  x12[T]. Also, we know that T = (d/(Subscript[k, F]))for this reason we make the replacement T -> (d/(Subscript[k, F]))

eq2 = 
        DSolve[{(Subscript[x, a])'[t] ==  - (Subscript[k, a]) (Subscript[x, a])[t], x1'[t] == 
              (Subscript[k, a]) (Subscript[x, a])[t] -  ((Subscript[k, 10]) + (Subscript[k, 12])) x1[t] + 
                (Subscript[k, 21]) x2[t], x2'[t] == (Subscript[k, 12]) x1[t] - (Subscript[k, 21]) x2[t], 
            (Subscript[x, a])[0] == (Subscript[X, a]), x1[0] == X1, 
            x2[0] == X2}, {(Subscript[x, a])[t], x1[t], x2[t]} ,  t];

{x2a[t_, d_], x21[t_, d_], 
          x22[t_, d_]}  = {(Subscript[x, a])[t], x1[t], x2[t]} /. 
              eq2[[1]]  /. {(Subscript[X, a]) ->  x1a[T], 
              X1 ->  x11[T], X2 ->  x12[T]} /. 
          T ->  (d/(Subscript[k, F])) ;

Now we will build a compound function that represent the retention for both intervals: one for  0 < t ≤ T where the input is constant and the other for t  > T when there isn't any input. We apply If whose syntax is If[condition, A, B] gives A if condition evaluates to True, and B if it evaluates to False. In our case A correspond to the solutions for the first interval given by eq13 (condition t ≤ T), and B the solutions for the second interval given by eq14 (condition t > T).

xa[t, d], x1[t,d]  and x2[t,d], represent the retentions in compartment a, 1 and 2. We live d as function variable better than T because T is a function of d.

xa[t_, d_] :=  
      If[ t ≤ (d/(Subscript[k, F])), x1a[t],  x2a[t - (d/(Subscript[k, F])), d]]

x1[t_, d_] :=  
      If[ t ≤ (d/(Subscript[k, F])), x11[t],  x21[t - (d/(Subscript[k, F])), d]]

x2[t_, d_] :=  
      If[ t ≤ (d/(Subscript[k, F])), x12[t],  x22[t - (d/(Subscript[k, F])), d]]

Here is plot the retention. We have suppose d = 10

Plot[{xa[t, 10], x1[t, 10], x2[t, 10]}, {t, 0, 200}, 
      PlotRange ->  All]

Clear[xa, x1, x2, model1, eq1,  eq2 ];

The lead is one of typical pollutant in some industrial environment. It can be intaken to the human body specially through inhalation. The intaken lead goes to blood, to the tissues and to the bone. Part of it is eliminated by urine excretion. We can represente the temporal evolution in the human body using a system of tree compartment: (1) blood, (2) tissues, (3) bones. The flow diagram is shown in the figure

model  = 
      CompartMatrix[
        3,  {{1, 2, k12}, {2, 1, k21}, {1, 3, k13}, {3, 1, k31}, {1, 0,  
            k10}, {2, 0, k20}}]

It can be represented by the system of diferential equation

eq1 = 
        ShowODE[model,  {0, 0, 0}, {b[t], 0, 0},  t,  x];

eq1 // TableForm // TraditionalForm

We will use the habitual agreement calling the lead amount in compartment 1, 2, y 3,  in the time t,  by  (Subscript[x, 1])(t), (Subscript[x, 2])(t) and (Subscript[x, 3])(t).  We call (Subscript[b, 1])the lead intake from outside of system in compartment 1.

In matrix notation  the previous equation is transformed in

x'(t) =A x + b(t), t ≥ 0

with

x'(t) =(({{(Subscript[x, 1])'}, {(Subscript[x, 2])'}, {(Subscript[x, 3])'}})); A=(({{(Subscript[a, 11]), (Subscript[a, 12]), (Subscript[a, 13])}, {(Subscript[a, 21]), (Subscript[a, 22]), 0}, {(Subscript[a, 31]), 0, (Subscript[a, 33])}})); x=(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}}));b(t)= (({{(Subscript[b, 1])(t)}, {0}, {0}}))

where the coefficiente (Subscript[a, ij]) of matrix A are:

(Subscript[a, 11]) =  -((Subscript[k, 10]) + (Subscript[k, 12]) + (Subscript[k, 13])),  
    (Subscript[a, 12]) = (Subscript[k, 21]),  (Subscript[a, 13]) = (Subscript[k, 31]),   (Subscript[a, 21]) = (Subscript[k, 12]),  
    (Subscript[a, 22]) =  -((Subscript[k, 20]) + (Subscript[k, 21])),  (Subscript[a, 31]) = (Subscript[k, 13]),  
    (Subscript[a, 33]) = -(Subscript[k, 31])

To resolve the previous system we must know the form of input b(t), also of the inital conditions x (0)= (Subscript[x, 0])

In the model of figure we consider a single intake (Subscript[b, 0]) in t  =  0 , so b(t) = (Subscript[b, 0]) for t  =  0 and b(t)  =  0 for t > 0. This is equivalent to assume as initial (Subscript[x, 1])(0) = (Subscript[b, 0]),  and (Subscript[(x ), 2])(0) = 0. If we try to use DSolve without giving values to (Subscript[a, ij]) we will see that no solution is obtain. In real problems the most usual thing is that transfer coefficient values are known.

In this example we will use the values given by 1992 (Bo92), in (days^(-1)), obtained empirically from a person exposure to an environment highlly contaminated intaken 49.3 μg of lead every day

(Subscript[k, 21])  =  0.0124; (Subscript[k, 12]) = 0.0111; (Subscript[k, 13]) = 0.0039; 
    (Subscript[k, 31])  =  0.000035; (Subscript[k, 10])  =  0.0211; 
    (Subscript[k, 20]) =  0.0162; 
Initial conditions :  b1 = 49.3   ,  b2  = 0;  
    b3 = 0;

model1 = 
        model /. {k21  ->  0.0124, k12 -> 0.0111, 
            k13 -> 0.0039,   k31  ->  0.000035,  k10  ->  0.0211,  
            k20 ->   0.0162};

{x1[t1_], x2[t1_], x3[t1_]} = {(Subscript[x, 1])[t1], (Subscript[x, 2])[t1], (Subscript[x, 3])[t1]} /. 
        SystemDSolve[model1, {0, 0, 0}, {49.3, 0, 0},  t,  t1,  
          x]

Plot[{ x1[t], x2[t], x3[t]},  {t,  0,  400}, 
 
        PlotStyle -> {
{AbsoluteThickness[0.5],  
              AbsoluteDashing[{4, 4}]}, {AbsoluteThickness[
                1]}, {AbsoluteThickness[1.5]}},  
        AxesLabel -> {"days",  "r"}, 
        PlotLegend -> {"blood", "tissues", "bones"}, 
        LegendPosition -> {1,  0}, 
        LegendSize -> {0.6,  
            0.5}, 
PlotLabel -> "Lead model:input = \
49.3 μg/day"];

Clear[b1, k, a, x1, x2, x3,  eq1,  model];

The figure represent the general tricompartmental model with one input in each compartment

This model is represented by the follow ODE System

(Subscript[x, 1])' (t) = (Subscript[b, 1])(
          t)  -  ((Subscript[k, 10])  +  (Subscript[k, 12])  +  (Subscript[k, 13])) (Subscript[x, 1])(
            t)   +  
        (Subscript[k, 21]) (Subscript[x, 2])(t)  +  (Subscript[k, 31]) (Subscript[x, 3])(
            t) 
(Subscript[x, 2])' (t) = (Subscript[b, 2])(t)  +  
        (Subscript[k, 12]) (Subscript[x, 1])(
            t)   -  ((Subscript[k, 20])  +  (Subscript[k, 21])  +  (Subscript[k, 23])) (Subscript[x, 2])(
            t)  +  (Subscript[k, 32]) (Subscript[x, 3])(t)
(Subscript[x, 3])' (t) = (Subscript[b, 3]) (t)  +  (Subscript[k, 13]) (Subscript[x, 1])(t)   +  
        (Subscript[k, 23]) (Subscript[x, 2])(
            t)  -  ((Subscript[k, 30])  +  (Subscript[k, 31])  +  (Subscript[k, 32])) (Subscript[x, 3])(
            t)

in matrix notation

x'(t) =A x + b(t), t ≥ 0

with

x\(TraditionalForm\`\((t)\) = \ {\(x\_1\)(t), \ \(x\_2\)(t), \ \(x\_1\)(t)}\^T\) 

B(t) = {(Subscript[b, 1])(t),  (Subscript[b, 2])(
                t), (Subscript[b, 3])(t)}^T 

where  A is

model  = 
        CompartMatrix[
          3,  {{1, 2, (Subscript[k, 12])}, {2, 1, (Subscript[k, 21])}, {1, 3, (Subscript[k, 13])}, {3, 1, 
              (Subscript[k, 31])}, {3, 2, (Subscript[k, 32])}, {2, 3, (Subscript[k, 23])}, {3, 3, (Subscript[k, 33])}, {1, 0,  
              (Subscript[k, 10])}, {2, 0, (Subscript[k, 20])}, {3, 0, (Subscript[k, 30])}}];

model // MatrixForm

It assumed as initial condition:

ic  = {(Subscript[x, 10]),  (Subscript[x, 20]),  (Subscript[x, 30])};

and as input

input = {(Subscript[b, 1]),  (Subscript[b, 2]),  (Subscript[b, 3])};

It can be represented by the system of diferential equation

eq1 = ShowODE[model,  ic, input,  t,  x];

eq1 // TableForm // TraditionalForm

We consider (Subscript[b, 1])(t)  =  2 Exp[-0.1 t], (Subscript[b, 2])(t) = 2,  (Subscript[b, 3])(t)  = 0.2  Cos[2 t] that is equivalent to take as initial conditions (Subscript[x, 1])(0)  = 1, (Subscript[x, 2])(0) = 2 , (Subscript[x, 3])(0)  =  3. Also we give values for constant coefficients.

k /: (Subscript[k, 10]) = 0.1; k /: (Subscript[k, 20]) = 0.2; k /: (Subscript[k, 30]) = 0.5; k /: (Subscript[k, 12]) = 1;  
    k /: (Subscript[k, 13]) = 2;  k /: (Subscript[k, 21]) = 1.4; k /: (Subscript[k, 23]) = 1.1;;  
    k /: (Subscript[k, 31]) = 1.1; k /: (Subscript[k, 32]) = 0.8;

Here is the solution

{x1[t1_], x2[t1_], x3[t1_]} = {(Subscript[x, 1])[t1], (Subscript[x, 2])[t1], (Subscript[x, 3])[t1]} /. 
        SystemDSolve[
          model, {1, 2, 3}, {2 Exp[-0.1 t], 2, 0.2  Cos[2 t] },  
          t,  t1,  x]

Plot[{x1[t], x2[t], x3[t]},  {t, 0, 100}, 
 
        PlotStyle -> {{AbsoluteThickness[2]}, {AbsoluteThickness[
                0.5]}, {AbsoluteThickness[1.5]}},  
        AxesLabel -> {"days",  "r"}, 
        PlotLegend -> {"x1", "x2", "x3"}, 
        LegendPosition -> {0,  0.3},  
        LegendSize -> {0.3,  0.3}];

Clear[model, k, input,  ic, x1, x2, x3, eq1];

Here is show the metabolic uranium model as appeared in ICRP 78

The compartments has been numbered as follow:

({{1, "Plasma"}, {2, "RBC"}, {3, "ST0"}, {4, "ST1"}, {5, "ST2"}, {6, "Liver1"}, {7, "Liver2"}, {8, "UrinaryPath"}, {9, "OtherKidney"}, {10, "TrabecularSurface"}, {11, "ExchTrabVol"}, {12, "NonExchTrabVol"}, {13, "CorticalSurf"}, {14, "CorticalVolExch"}, {15, "CorticalVolNonExch"}, {16, "Bladder"}, {17, "Urine"}, {18, "ULI"}, {19, "LLI"}, {20, "FEC"}})

Transfer rates are taken from Table A.12.1. (ICRP 78) (For GI compartments are applied the trasfer rate included in ICRP 30]).

({{"Plasma to ST0", "k[1,3]->", 10.5`}, {"Plasma to RBC", "k[1,2]->", 0.245`}, {"Plasma to urinary bladder", "k[1,16]->", 15.43`}, {"Plasma to Kidneys (urinary path)", "k[1,8]->", 2.94`}, {"Plasma to other kidney tissue", "k[1,9]->", 0.0122`}, {"Plasma to ULI contents", "k[1,18]->", 0.122`}, {"Plasma to liver 1", "k[1,6]->", 0.367`}, {"Plasma to ST1", "k[1,4]->", 1.63`}, {"Plasma to ST2", "k[1,5]->", 0.0735`}, {"Plasma to trabecular surface", "k[1,10]->", 2.04`}, {"Plasma to cortical surfaces", "k[1,13]->", 1.63`}, {"ST0 to plasma", "k[3,1]->", 8.32`}, {"RBC to Plasma", "k[2,1]->", 0.347`}, {"Kidneys (Urinary path) to urinary bladder", "k[8,16]->", 0.099`}, {"Other kidney tissue to plasma", "k[9,1]->", 0.00038`}, {"Liver 1 to plasma", "k[6,1]->", 0.092`}, {"Liver 1 to liver2", "k[6,7]->", 0.00693`}, {"ST1 to plasma", "k[4,1]->", 0.0347`}, {"ST2 to plasma", "k[5,1]->", 0.000019`}, {"Cortical(Bone) surfaces to plasma", "k[13,1]->", 0.0693`}, {"Cortical (Bone) surface to exch cortical volume", "k[13,14]->", 0.0693`}, {"Liver 2 to plasma", "k[7,1]->", 0.00019`}, {"Nonexch trab vol to plasma", "k[12,1]->", 0.000493`}, {"Nonexch cort vol to plasma", "k[15,1]->", 0.0000821`}, {"Cortical(Exch bone) vol to cortical (bone) surfaces", "k[14,13]->", 0.0173`}, {"Exch trabecular vol to nonexch trabecular vol", "k[11,12]->", 0.00578`}, {"Trabecular surface to plasma", "k[10,1]->", 0.0693`}, {"Exch cortical vol to nonexch cortical vol", "k[14,15]->", 0.00578`}, {"Exch trabecular vol to trabecular surface", "k[11,10]->", 0.0173`}, {"Trabecular surface to Exchan. Trabecular vol.", "k[10,11]->", 0.0693`}})

({{"Bladder (16) to urine(17)", "k[16,17]->", 12}, {"ULI (18) to LLI(19)", "k[18,19]->", 1.8}, {"LLI (19) to FEC(20)", "k[19,20]->", 1}})

The transfer coefficients for uranium are the following

uraniumModel = {{1, 3, 10.5`}, {1, 2, 
            0.245`}, {1, 16, 15.43`}, {1, 8, 2.94`}, {1, 9, 0.0122`}, {1, 18, 
            0.122`}, {1, 6, 0.367`}, {1, 4, 1.63`}, {1, 5, 0.0735`}, {1, 10, 
            2.04`}, {1, 13, 1.63`}, {3, 1, 8.32`}, {2, 1, 0.347`}, {8, 16, 
            0.099`}, {9, 1, 0.00038`}, {6, 1, 0.092`}, {6, 7, 0.00693`}, {4, 
            1, 0.0347`}, {5, 1, 0.000019`}, {13, 1, 0.0693`}, {13, 14, 
            0.0693`}, {7, 1, 0.00019`}, {11, 12, 0.00578`}, {12, 1, 
            0.000493`}, {15, 1, 0.0000821`}, {14, 13, 0.0173`}, {14, 15, 
            0.00578`}, {10, 1, 0.0693`}, {11, 10, 0.0173`}, {10, 11, 
            0.0693`}, {16, 17, 12}, {18, 19, 1.8}, {19, 20, 
            1}};

The instant urine and faecal rate excrection for a single input 1 in compartment 1 can be evaluated as follow

{qAcuteUrine[t_],  qAcuteFaecal[t_]} = 
      D[{(Subscript[x, 17])[t],  (Subscript[x, 20])[t]} /. 
          AcuteInput[CompartMatrix[20, uraniumModel], PadRight[{1}, 20] ,  
            t,   x], t]

LogPlot[{qAcuteUrine[t],  qAcuteFaecal[t]},  {t,  2, 10}]

The same problem asuming a constant intake (chronic) an single input 1 in compartment 1 can be evaluated as follow

{qConstantUrine[t1_],  qConstantFaecal[t1_]} = 
      D[{(Subscript[x, 17])[t1],  (Subscript[x, 20])[t1]} /. 
          SystemDSolve[CompartMatrix[20, uraniumModel], PadRight[{0}, 20] , 
            PadRight[{1}, 20],  t,  t1,  x],  t1]

LogPlot[{qConstantUrine[t],  qConstantFaecal[t]},  {t,  2, 
        100}]

Clear[uraniumModel, qAcuteUrine,  qAcuteFaecal,  qConstantUrine,  
        qConstantFaecal];

In this model is supposed that the transfer  (Subscript[k, bilis, gi]) from Bilis to GI it is a UnitStep function [h[t,8, 2]]

h[t_,  a_,  d_] :=1 - UnitStep[Mod[t,  a + d] - d];

Plot[h[t, 8,  2], {t, 0, 110}, PlotPoints -> 1000, 
        PlotRange -> All];

The clearance coeffs are{k12 -> 2, k21 ->  1, k10 ->  0.5, 
        (Subscript[k, g]) ->  0.1, (Subscript[k, a]) ->  2, 
        (Subscript[k, bilis, gi]) ->  h[t, 8, 2] }

matrixhepatic = 
        CompartMatrix[
            4,  {{1, 2,  k12},  {1, 0, k10}, {2, 1,  k21}, {1, 3,  
                (Subscript[k, g])}, { 3, 4,  (Subscript[k, bilis, gi])},  {4, 1,  
                (Subscript[k, a])}}] /. {k12 -> 2, k21 ->  1, 
            k10 ->  0.5, (Subscript[k, g]) ->  0.1, (Subscript[k, a]) ->  2, 
            (Subscript[k, bilis, gi]) ->  h[t, 8, 2] };

Initial condition and input are

ip  = {0, 0, 0, 0};

ic2  = {200, 0, 0, 0};

Then, it can be used  SystemNDSolve

sol1 = 
      SystemNDSolve[matrixhepatic,  ic2, ip, { t, 0,  100},  
        t,   x]

We need increase the number of steps

sol2 =  
      SystemNDSolve[matrixhepatic,  ic2, ip, { t, 0,  100},  t,   x, 
        MaxSteps -> 5000]

{ x1[t_], x2[t_], x3[t_], 
        x4[t_]}  =  { (Subscript[x, 1])[t], (Subscript[x, 2])[t], (Subscript[x, 3])[t], (Subscript[x, 4])[t]} /. 
        sol2

Here are ploted the solution

Plot[{x1[t], x2[t], x3[t]},  {t, 0, 50}, PlotPoints -> 1000, 
        PlotRange -> {0, 100}, 
 
        PlotStyle -> {{AbsoluteThickness[2]}, {AbsoluteThickness[
                0.5]}, {AbsoluteThickness[1.5]}},  
        AxesLabel -> {"days",  "r"}, 
        PlotLegend -> {"x1", "x2", "x3"}, 
        LegendPosition -> {0,  0.3},  
        LegendSize -> {0.3,  0.3}];

Here are shown the numerical values

Table[{t,  x1[t],  x2[t],  x3[t],  x4[t]},  {t,  1,  100,  5}] // 
      TableForm

The one compartmental  Michaeles-Menten model is given by equation

\(x\& . \_1\)(t) = - (((Subscript[V, m]) (Subscript[x, 1])[t])/((Subscript[k, m]) + (Subscript[x, 1])[t]))+ (Subscript[b, 1]) (Subscript[u, 1])[t]

(Subscript[b, 1]) (Subscript[u, 1])[t]= (Subscript[D, 1]) δ[t], this men that the initial condion is   (Subscript[x, 1])[0] =  (Subscript[D, 1])

DSolve[{(Subscript[x, 1])'[t] + (((Subscript[V, m]) (Subscript[x, 1])[t])/((Subscript[k, m]) + (Subscript[x, 1])[t])) == 
          0,  (Subscript[x, 1])[0]  ==  (Subscript[D, 1])}, (Subscript[x, 1])[t],  t]

% // FullSimplify

m[t_] =  (Subscript[x, 1])[t] /. %

Plot[
      m[t] /. {(Subscript[D, 1]) ->  1, (Subscript[k, m])  ->  12,  
          (Subscript[V, m])  ->  15},  {t,  0, 2}]

Clear[x1, x2]

The bicompartmental  Michaeles-Menten model is given by equation

\(x\& . \_1\)(t) =  - (Subscript[k, 21]) (Subscript[x, 1])(t)  + (Subscript[k, 12]) (Subscript[x, 2])(t) -  \(TraditionalForm\`\(\(\ \)\(\(k\_1\)\(\ \)\(\(x\_1\)(t)\)\(\ \)\)\)\/\(k\_2 + \ \(x\_1\)(t)\)\) 
\(x\& . \_2\)(t) = (Subscript[k, 21]) (Subscript[x, 1])(t)  -  (Subscript[k, 12]) (Subscript[x, 2])(t) 
(Subscript[x, 1])(0) = (Subscript[b, 0]) , (Subscript[x, 2])(0) = 0

with clearance coeffs:

(Subscript[k, 1])  =  0.2;  (Subscript[k, 2])  =  0.3; (Subscript[k, 12])  =  0.03; 
    (Subscript[k, 21])  =  0.02;

eq1  =  x1'[t] + (Subscript[k, 21]) x1[t] - 
          (Subscript[k, 12]) x2[t] + (((Subscript[k, 1]) x1[t])/((Subscript[k, 2]) + x1[t])) == 0

eq2  =  x2'[t] - (Subscript[k, 21]) x1[t] + (Subscript[k, 12]) x2[t] == 0

It can be solved using NDSolve

NDSolve[{eq1,  eq2, x1[0]  ==  1, x2[0]  ==  0}, {x1[t], 
        x2[t]},  {t, 0, 10}]

Plot[Evaluate[{x1[t], x2[t]} /. %], {t, 0, 10}, 
        PlotPoints -> 100, 
        PlotRange -> 
          All,   PlotStyle -> {{AbsoluteThickness[
                2]}, {AbsoluteThickness[0.5]}}, 
        PlotLegend -> {"x1", "x2"}, 
        LegendPosition -> {0,  0.2}];

k  =.

Clear[eq1, eq2]

Here is solved eq. (1) by Roberts-Rowland J. Pharmacokim.Biopharm 14: 227-260 (1986)

\(∂\^2 c\/∂z\^2 - \(1\/d\_n\) ∂c\/∂z - \(r\_n\/d\_n\) c = \(1\/d\_n\) ∂c\/∂z, 
    where\ c = \(\(c\) \((z, t)\)\(\ \)\)\)

Here is solved the case with "mixed" boundary conditions

Initial condition (IC): t =0, c= 0, 0≤z≤1
Boundary conditions (BC): (The "mixed" boundary conditions): For z=0: c =δ(z), for z=∞: c=0, in other notation: c(0) = δ(t) y c(∞)  = 0.

Eq(17) in Mathematica notation can be wrote

D[c[z, t], t] == 
      (Subscript[d, n]) D[c[z, t], z, z] -  D[c[z, t], z] - 
        (Subscript[r, n]) (Subscript[d, n]) c[z, t] (2)

IC: c[z,0]->0

Here is solved eq.(18) using Laplace Transform. (The right side terms of eq. 18 are denoted by eq1 and the left by eq2)

eq1 = LaplaceTransform[D[c[z, t], t],  t, s]

eq2 = 
      LaplaceTransform[
        (Subscript[d, n]) D[c[z, t], z, z] -  D[c[z, t], z] - (Subscript[r, n]) (Subscript[d, n]) c[z, t],  t, 
        s]

Now we calculate their Laplace Transforms. The are denoted in capital letters (i.e: for c[z,0] the LT is C[z,s] ). BC and IC are applied.

eq1LT  =  
      eq1 /. {LaplaceTransform[c[z, t], t, s] -> C[z, s]}

eq2LT = eq2/.{LaplaceTransform[Derivative[1, 0][c][z, t],t,s]->Derivative[1, 0][C][z, s],LaplaceTransform[Derivative[2, 0][c][z, t],t,s]->Derivative[2, 0][C][z, s],LaplaceTransform[c[z, t], t, s] -> C[z, s]}

A26 = eq1LT  == eq2LT /. {c[z, 0] -> 0}

LaplaceTransform[DiracDelta[z], z, s]

eq26b = 
      DSolve[{s C[z] == -C[z] (Subscript[d, n]) (Subscript[r, n]) - C'[z] + 
              (Subscript[d, n]) C''[z], C[0] == 1}, C[z], z]

The integration constant C[2] can be calculated using the BC: (Subscript[C, b])(∞)  = 0

Solve[eq26b[[1]][[1, 2]] == 0,  
        C[2]];

c = C[2] /. %[[1]] // Simplify

eq26b /. {C[2] -> 0} // FullSimplify

Cout[s_] = C[z] /. %[[1]]

Now it can be evaluated  (Subscript[C, out]) (z, t)  taken the inverse laplace transform

cOut[z_, t_] = 
      InverseLaplaceTransform[Cout[s], s, t] // FullSimplify

When  z=1 the solution is:

cOut[1, t]

This solution has a little difference with  eq.(49) Roberts-Rowland J. Pharmacokim.Biopharm 14: 227-260 (1986): (Subscript[r, n]) was removed of the exponent. This mistake has not consequences when (Subscript[r, n])=0. as happen in  Fig 3.(B) of Roberts-Rowland paper

Here is shown this Fig 3.(B) obtain with our method {(Subscript[d, n]) -> {0.01, 0.1, 0.2, 1, 2}, (Subscript[r, n]) -> 0

Plot[
      Evaluate[cOut[1, t] /. {(Subscript[d, n]) -> {0.01, 0.1, 0.2, 1, 2}, 
            (Subscript[r, n]) -> 0}], {t, 0, 4},  
      PlotRange -> {0, 4}]

Here is the 3-D representation as function of z and t

Plot3D[
      Evaluate[cOut[z, t] /. {(Subscript[d, n]) -> 0.01, (Subscript[r, n]) -> 0}], {z, 1, 
        2}, {t, 0.1, 4},  PlotRange -> {0, 4}]

Plot3D[
      Evaluate[cOut[z, t] /. {(Subscript[d, n]) -> 0.1, (Subscript[r, n]) -> 0}], {z, 0.1, 
        2}, {t, 0.1, 4},  PlotRange -> {0, 4}]

Clear[c]

Here is shown that the solution obtained with Mathematica is right

\(\(c[z_, 
          t_] = \(E\^\(\(-\(\((t - z)\)\^2\/\(4\ t\ d\_n\)\)\) - t\ d\_n\ r\_n\)\ z\)\/\(2\ \@π\ \@\(t\^3\/d\_n\%3\)\ d\_n\%2\);\)\)

D[c[z, t], t] == 
        (Subscript[d, n]) D[c[z, t], z, z] -  D[c[z, t], z] - (Subscript[r, n]) (Subscript[d, n]) c[z, t] // 
      Simplify

[1]  J. A. Jacquez,  Compartmental Analysis in Biology and Medicine, The University of Michigan Press, 1985.

[2] Anderson, D.H. Compartmental analysis and tracer kinetics. Berlin: Springer-Verlag; 1983.

[3] Godfrey, K. Compartmental models and their application. London: Academic Press; 1983.