Sim

Initialization

Load the MCA package

<< MetabolicControlAnalysis`

Volumes

The total sample volume of the incubations was 2 ml and the cell volume was 19.3 microL per mL of incubation medium:

sampleV = 2.0*10^-3;
cellV = 38.68*10^-6;

The mitochondrial volume was 28% of the cell volume (10.8304 μL):

mitoV = 0.28*cellV;

The nucleus volume was 8% of the cell volume:

nucV = 0.08*cellV;

Thus the cytoplamic volume was:

cytV = cellV - mitoV - nucV;

The ratio of the imsV to total mitochondrial volume was 1.82/(1.82+2.99)= 37.8%. Hence the matrix volume was 62.2% of the mitochondrial volume (note that ims and mat denote 'inter-membrane space' and 'matrix', respectively):

imsV =  0.378*mitoV;
matV =  0.622*mitoV;

The volume of the 'extra mitochondrial' space is therefore:

eV = sampleV - cellV;

General Initial Conditions

These are the total enzyme concentrations expressed per kg (or L) of liver. But now changed to give the corresponding concentration in the cytoplasm or mitochondial matrix, as required. This is according to PWK's original urea cycle simulation (1977):

e10 := (cellV/
              matV)*2.6 10^-6;  
e20 := (cellV/
              cytV)*4.0 10^-6;  
e30 := (cellV/
              cytV)*2.2 10^-6;

(*calculate the concentration of arginase inside and outside the cells*) 

e40origConc = 
        8.9*10^-6; (*conc of arginase mol/
          kg liver*) 


    e40totMol = 
        e40origConc*
          cellV; (*amount of arginase in a cell (mol)*) 


    arginaseLeakage = 
        0.13; (*proportion of arginase leaked out of the cell, 
      according to the methods it was about 13  %*) 


    e40inAmt = (1 - arginaseLeakage)*
          e40totMol; (*amount of arginase in the cytoplasm*) 


    e40outAmt = 
        arginaseLeakage*
          e40totMol; (*amount of arginase in the extracellular medium*) 


    e40 = 
      e40inAmt/cytV (*concentration of arginase in the cytoplasm*) 


    e40E = e40outAmt/
        eV (*concentration of arginase in the extracellular medium*)

These are the input pool concentrations:

atppool = 2.0 10^-3;
asppool = 1.0 10^-3;
cppool = 15*10^-3;

General Initial Conditions

The 'specific activity': (NOTE: here we use carbamoyl phosphate as the added radioactive species.)
Here we recall that in the experiments, the total amount of bicarbonate was 15 mM. This translates to 30 μmol in the 2 ml suspensions. However, the ammonia concentration was 0.5 mM in the experiment we are concerned with. Assuming this rapidly equilibrates across the mitochondrial membrane, it's total concentration will also be 0.5 mM in the matrix. The means that the total potential concentration of the cp in the matrix can only be as high as 0.5 mM. So we will assume that this is the initial, total concentration of cp in the matrix. We will also assume that the radiolabelled bicarbonate bought by the authors was 100% labelled, that is, every molecule contained 14C. According to the methods section, this had a specific activity of 55 mCi/mmol. This is equal to 1.221*10^8 dpm/μmol. Now if we know that the specific activity of the bicarbonate in experiment 1 was 530 cpm/nmol (= 15900000 cpm/30 μmol), and knowing that the counting efficiency for 14C was 96%, we can calculate that there must have been 0.136 μmol of radiolabelled bicarbonate in every 30 μmol of this metabolite. In other words every mole of carbamoyl phosphate would have 0.00453333 moles of radioactive molecules. Hence we now have a starting point for the initial concentration of labelled intermediates.

sa = 530;

ratio = 
        0.004533333;

cpZ := 
      0.5*10^-3;  (*initial carbamol phosphate concentration in the extracellular space*) 
cpRZ := cpZ*ratio (*initial RADIOACTIVE carbamol phosphate concentration in the extracellular space*)

Input the initial concentrations of each metabolite:

initialconditions = {
cp[0. ] == 
            1*10^-5, 
cpe[0. ] == 
            1*10^-5, 
cpeE[0] == 
            cpZ, 
cpims[0. ] == 
            1*10^-5, 
cpeR[0. ] == 
            0, 
cpeRE[0] == 
            cpRZ, 
cpimsR[0. ] == 
            0, 
cpR[0. ] == 
            0, 
ornmat[0. ] == 
            1*10^-6, 
citmat[0. ] == 
            1*10^-6, 
citmatR[0. ] == 
            0, 
pi[0. ] == 
            1.0 10^-4, 
citims[0. ] == 
            1*10^-6, 
citimsR[0. ] == 
            0, 
cit[0. ] == 
            1*10^-6, 
citR[0. ] == 
            0, 
asp[0. ] == 
            1.0*10^-4, 
atp[0. ] == 
            1.5*10^-3, 
as[0. ] == 
            1*10^-6, 
asR[0. ] == 
            0, 
amp[0. ] == 
            1*10^-6, 
ppi[0. ] == 
            1*10^-6, 
fum[0. ] == 
            1*10^-6, 
arg[0. ] == 
            1*10^-6, 
argQ[0. ] == 
            0, 
argR[0. ] == 
            0, 
orn[0. ] == 
            0.2*10^-3, 
ornims[0. ] == 
            4.5 10^-4, 
ure[0.0] == 
            1*10^-6, 
ureR[0.0] == 
            0, 
citE[0] == 
            1*10^-6, 
citRE[0] == 
            0, 
atpE[0] == 
            1.5*10^-3, 
aspE[0] == 
            1*10^-4, 
ampE[0] == 
            1*10^-6, 
ppiE[0] == 
            1*10^-6, 
asE[0] == 
            1*10^-6, 
asRE[0] == 
            0, 
fumE[0] == 
            1*10^-6, 
argQE[0] == 
            argAdd, 
argE[0] == 
            1*10^-6, 
argRE[0] == 
            0, 
ureE[0] == 
            1*10^-6, 
ureRE[0] == 
            0, 
ornE[0] == 
            0.2*10^-3
};

(*A note on the naming of metabolites:If the metabolite name is suffixed with 'ims', 'mat' or 'E' this denotes that it is in the inter-membrane space, mitochondrial matrix,or the extracellular space,respectively.Otherwise it is generally belonging to the cytoplasm.If it is suffixed with or contains an upper-case 'R' it is RADIOACTIVE.
 Exceptions are matrix carbamoyl phosphate (which is just cp[t]) and matrix inorganic phosphate (pi[t])*)

Unitary Rate Constants

For information on the derivation of these see appendix 1.

Ornithine carbamoyltransferase:

k1e1 = 1.7 10^7;
k2e1 = 6.3 10^1;
k3e1 = 2.1 10^6;
k4e1 = 1.0 10^3;
k5e1 = 3.0 10^3;
k6e1 = 9.0 10^4;
k7e1 = 2.6 10^3;
k8e1 = 5.0 10^5;

Argininosuccinate synthase:

k1e2 = 2.4 10^5;
k2e2 = 2.3;
k3e2 = 3.5 10^5;
k4e2 = 10.0;
k5e2 = 4.8 10^5;
k6e2 = 10.0;
k7e2 = 2.0 10^1;
k8e2 = 8.9 10^5;
k9e2 = 5.0 10^1;
k10e2 = 6.4 10^5;
k11e2 = 5.0 10^1;
k12e2 = 1.7 10^5;

Argininosuccinate lyase:

k1e3 = 2.7 10^6;
k2e3 = 7.0 10^1;
k3e3 = 7.5 10^1;
k4e3 = 1.5 10^6;
k5e3 = 1.1 10^3;
k6e3 = 7.0 10^5;

Arginase:

k1e4 = 1.0 10^7;
k2e4 = 5.4 10^4;
k3e4 = 5.3 10^3;
k4e4 = 3.0 10^4;
k5e4 = 1.0 10^7;

Rate constants for input or efflux of 'external' metabolites:

\(\(\( (*the\ first\ step\ is\ to\ get\ a\ number\ for\ the\ rate\ of\ production\ of\ urea\ by\ the\ cells\ if\ there\ was\ a\ ' steady\ state' 
        in\ units\ of\ mol/s*) \)\(
\)\(\(uppm = 
        12.1;\)\(\ \)\( (*the\ rate\ of\ production\ of\ urea\ in\ the\ experiments, \(nmol/min\)/ml\ of\ cells*) \)\(
\)
    \(uppm1 = uppm*10^\(-9\);\)\( (*\(mol/mL\)/min*) \)\(
\)
    \(uppm2 = uppm1*1000;\)\( (*\(mol/L\)/min*) \)\(
\)
    \(uppm3 = uppm2/60;\)\(\ \)\( (*\(mol/L\)/s*) \)\(
\)
    \(uppm4 := uppm3*cellV;\)\(\ \)\( (*mol/s*) \)\(
\)
    \(steadyStateUreaProdn := uppm4;\)\(\ \)\( (*mol/
        s*) \)\(
\)\(
\)\( (*now\ we\ can\ define\ the\ rate\ constants\ for\ input\ of\ the\ peripheral\ metabolites\ into\ their\ respective\ compartments\ to\ be\ identical\ \((in\ mol/s)\)\ to\ the\ ' 
        steady\ state'\ rate\ of\ urea\ production\ defined\ above*) \)\(
\)
    \(katp := 
      steadyStateUreaProdn/\((cytV*
            atppool)\);  (*rate\ constant\ for\ ATP\ input*) 

    kasp := steadyStateUreaProdn/\((cytV*
            asppool)\);\)\( (*rate\ constant\ for\ aspartate\ input*) \)\(
\)
    \(kcp := 
        steadyStateUreaProdn/\((eV*
              cppool)\);\)\( (*rate\ constant\ for\ carbamoyl\ phosphate\ input*) \)\(
\)\(
\)\( (*these\ are\ the\ desired\ final\ concentrations\ for\ the\ peripheral\ metabolites\ produced\ in\ the\ cycle
All\ in\ mol/L*) \)\(
\)
    \(ppifinal = 
        0.01*10^\(-3\);\)\(\ \)\( (*steady\ state\ pyrophosphate\ concentration*) \)\(
\)
    \(fumfinal = 
        0.01*10^\(-3\);\)\( (*steady\ state\ fumarate\ concentration*) \)\(
\)
    \(ampfinal = 
        0.05*10^\(-3\);\)\( (*steady\ state\ AMP\ concentration*) \)\(
\)
    \(pifinal = 
        0.05*10^\(-3\);\)\( (*steady\ state\ inorganic\ phosphate\ in\ the\ matrix\ concentration*) \)\(
\)
    \(urefinal = 
        0.001*10^\(-3\);\)\( (*steady\ state\ urea\ concentration*) \)\(
\)\(
\)\( (*and\ we\ can\ now\ define\ the\ rate\ constants\ for\ REMOVAL\ of\ peripheral\ metabolites\ from\ their\ respective\ compartments\ to\ be\ identical\ to\ the\ desired\ ' 
        steady\ state'\ rate\ of\ production\ of\ urea*) \)\(
\)
    \(kppi := 
        steadyStateUreaProdn\ /\((cytV*
              ppifinal)\);\)\(\ \)\( (*rate\ constant\ for\ pyrophosphate\ removal*) \)\(
\)
    \(kfum := 
        steadyStateUreaProdn/\((cytV*
              fumfinal)\);\)\( (*rate\ constant\ for\ fumarate\ removal*) \)\(
\)
    \(kamp := 
        steadyStateUreaProdn/\((cytV*
              ampfinal)\);\)\( (*rate\ constant\ for\ AMP\ removal*) \)\(
\)
    \(kpi := 
        steadyStateUreaProdn/\((matV*
              pifinal)\);\)\( (*rate\ constant\ for\ inorganic\ phosphate\ removal*) \)\(
\)
    \(kure := 
        steadyStateUreaProdn/\((cytV*
              urefinal)\);\)\( (*rate\ constant\ for\ urea\ removal*) \)\(
\)\(
\)
    \( (*note\ that\ the\ above\ are\ input\ to/
          removed\ from\ the\ CYTOPLASM\ even\ though\ they\ are\ freely\ permeable\ to\ the\ extracellular\ medium . \ \((From\ a\ simulation\ point\ of\ view\ it\ would\ make\ no\ difference\ if\ they\ were\ input\ to/
                removed\ from\ the\ extracellular\ space, 
            for\ example)\)*) \)\)\)\)

Ki of arginine for inhibition of arginosuccinate synthase (see 3rd paragraph of column 2 p 4040 of the Cheung JBC paper)

Kiarg = 0.7 10^-3;

Rate constants for membrane exchange reactions,viz.,transporters.

kornci := 1*10^4;

kornic := kornci cytV/imsV;

kcitci = 1*10^4;

kcitic := kcitci cytV/imsV;

kornim = 1.1*10^2;

kornmi := kornim imsV/matV;

kcitim = 1*10^2;

kcitmi := kcitim imsV/matV;

kcpei = 1 10^4;

kcpie := kcpei cytV/imsV;

kcpmi = 0.5 10^0;

kcpim := kcpmi matV/imsV;

(*the following are the unitary rate constants for transport of metabolites from the cytoplasm to the extracellular medium . All are given a high value to reflect fast transport*) 

kcitCE := 
        1*10^7;

    kcitRCE  := 1*10^7;

    katpCE := 1*10^7;

    kaspCE := 1*10^7;

    kampCE := 1*10^7;

    kppiCE := 1*10^7;

    kasCE := 1*10^7;

    kasRCE := 1*10^7;

    kfumCE := 1*10^7;

    kargQCE := 1*10^7;

    kargCE := 1*10^7;

    kargRCE := 1*10^7;

    kureCE := 1*10^7;

    kureRCE := 1*10^7;

    kornCE := 1*10^7;

    kcpeCE := 1*10^7;

    kcpeRCE := 1*10^7;

kcitCEREV := 
      1*10^7*cytV/eV

    kcitRCEREV := 1*10^7*cytV/eV

    katpCEREV := 1*10^7*cytV/eV

    kaspCEREV := 1*10^7*cytV/eV

    kampCEREV := 1*10^7*cytV/eV

    kppiCEREV := 1*10^7*cytV/eV

    kasCEREV := 1*10^7*cytV/eV

    kasRCEREV := 1*10^7*cytV/eV

    kfumCEREV := 1*10^7*cytV/eV

    kargQCEREV := 1*10^7*cytV/eV

    kargCEREV := 1*10^7*cytV/eV

    kargRCEREV := 1*10^7*cytV/eV

    kureCEREV := 1*10^7*cytV/eV

    kureRCEREV := 1*10^7*cytV/eV

    kornCEREV := 1*10^7*cytV/eV

    kcpeCEREV := 1*10^7*cytV/eV

    kcpeRCEREV := 1*10^7*cytV/eV

Rate Equations

Ornithine carbamoyltransferase:

otcAdj = otcAdj;

vel1FOR := 
          otcAdj*(k1e1 k3e1 k5e1 k7e1 cp[t] ornmat[t])*
            e10/denom1;

vel1REV := 
          otcAdj*(k2e1 k4e1 k6e1 k8e1 citmat[t] pi[t])*
            e10/denom1;

(*Velocity of the FORWARD ornithine carbamoyltransferase reaction when the cp IS radioactive.*) 
vel1FORRc := 
          otcAdj*(k1e1 k3e1 k5e1 k7e1 cpR[t] ornmat[t])*
            e10/denom1;

vel1REVRc := 
          otcAdj*(k2e1 k4e1 k6e1 k8e1 citmatR[t] pi[t])*
            e10/denom1;

denom1 := (k2e1 k7e1 (k4e1 + k5e1) + 
            k1e1 k7e1 (k4e1 + k5e1) (cp[t] + cpR[t]) + 
            k2e1 k8e1 (k4e1 + k5e1) pi[t] + 
            k3e1 k5e1 k7e1 ornmat[t] + 
            k2e1 k4e1 k6e1 (citmat[t] + citmatR[t]) + 
            k1e1 k3e1 (k5e1 + k7e1) (cp[t] + cpR[t]) ornmat[t] + 
            k6e1 k8e1 (k2e1 + k4e1) pi[
                t] (citmat[t] + citmatR[t]) + 
            k1e1 k4e1 k6e1 (cp[t] + cpR[t]) (citmat[t] + 
                  citmatR[t]) + 
            k1e1 k3e1 k6e1 (cp[t] + cpR[t]) 
              ornmat[t] (citmat[t] + citmatR[t]) + 
            k3e1 k5e1 k8e1 (ornmat[t]) pi[t] + 
            k3e1 k6e1 k8e1 (ornmat[t]) pi[
                t] (citmat[t] + citmatR[t]));

Argininosuccinate Synthase:

vel2FOR := (k1e2 k3e2 k5e2 k7e2 k9e2 k11e2 cit[t] atp[t] asp[t])*assAdj*
            e20*(1 + (arg[t] + argR[t] + argQ[t])/Kiarg)^-1/
              denom2;

(*Velocity of the REVERSE argininosuccinate synthase reaction when the argininosuccinate is NOT radioactive.*) 
vel2REV := (k2e2 k4e2 k6e2 k8e2 k10e2 k12e2 ppi[t] amp[
                  t] as[t])*assAdj*
            e20*(1 + (arg[t] + argR[t] + argQ[t])/Kiarg)^-1/
              denom2;

vel2FORRc := (k1e2 k3e2 k5e2 k7e2 k9e2 k11e2 citR[t] atp[t]*asp[t])*assAdj*
            e20*(1 + (arg[t] + argR[t] + argQ[t])/Kiarg)^-1/
              denom2;

vel2REVRa := (k2e2 k4e2 k6e2 k8e2 k10e2 k12e2 ppi[
                  t] amp[t] asR[t])*assAdj*
            e20*(1 + (arg[t] + argR[t] + argQ[t])/Kiarg)^-1/
              denom2;

denom2 := (k2e2 k4e2 k9e2 k11e2 (k6e2 + k7e2) + 
            k1e2 k4e2 k6e2 k8e2 k11e2 (cit[t] + citR[t]) ppi[t] + 
            k1e2 k4e2 k9e2 k11e2 (k6e2 + k7e2) (cit[t] + 
                  citR[t]) + 
            k2e2 k5e2 k7e2 k9e2 k12e2 asp[t] (as[t] + asR[t]) + 
            k2e2 k5e2 k7e2 k9e2 k11e2 asp[t] + 
            k1e2 k3e2 k6e2 k8e2 k11e2 (cit[t] + citR[t]) atp[
                t] ppi[t] + 
            k1e2 k3e2 k9e2 k11e2 (k6e2 + k7e2) (cit[t] + 
                  citR[t]) atp[t] + 
            k1e2 k4e2 k6e2 k8e2 k10e2 (cit[t] + citR[t]) ppi[
                t] amp[t] + 
            k1e2 k5e2 k7e2 k9e2 k11e2 (cit[t] + citR[t]) asp[t] + 
            k3e2 k5e2 k7e2 k9e2 k12e2 atp[t] asp[
                t] (as[t] + asR[t]) + 
            k3e2 k5e2 k7e2 k9e2 k11e2 atp[t] asp[t] + 
            k2e2 k5e2 k7e2 k10e2 k12e2 asp[t] amp[
                t] (as[t] + asR[t]) + 
            k1e2 k3e2 k5e2 (k7e2 k9e2 + k7e2 k11e2 + 
                  k9e2 k11e2)*(cit[t] + citR[t]) atp[t] asp[t] + 
            k1e2 k3e2 k5e2 k8e2 k11e2 (cit[t] + citR[t]) atp[
                t] asp[t] ppi[t] + k2e2 k4e2 k6e2 k8e2 k11e2 ppi[t] + 
            k1e2 k3e2 k5e2 k7e2 k10e2 (cit[t] + citR[t]) atp[
                t] asp[t] amp[t] + 
            k2e2 k4e2 k9e2 k12e2 (k6e2 + k7e2) (as[t] + 
                  asR[t]) + 
            k1e2 k3e2 k6e2 k8e2 k10e2 (cit[t] + citR[t]) atp[
                t] ppi[t] amp[t] + 
            k2e2 k4e2 k6e2 k8e2 k10e2 ppi[t] amp[t] + 
            k3e2 k5e2 k7e2 k10e2 k12e2 atp[t] asp[t] amp[
                t] (as[t] + asR[t]) + 
            k2e2 k4e2 k6e2 k8e2 k12e2 ppi[t] (as[t] + asR[t]) + 
            k3e2 k6e2 k8e2 k10e2 k12e2 atp[t] ppi[t] amp[
                t] (as[t] + asR[t]) + 
            k2e2 k4e2 k10e2 k12e2 (k6e2 + k7e2) amp[
                t] (as[t] + asR[t]) + 
            k2e2 k5e2 k8e2 k10e2 k12e2 asp[t] ppi[t] amp[
                t] (as[t] + asR[t]) + 
            k8e2 k10e2 k12e2 (k2e2 k4e2 + k2e2 k6e2 + 
                  k4e2 k6e2) ppi[t] amp[t]*(as[t] + asR[t]) + 
            k1e2 k3e2 k5e2 k8e2 k10e2 (cit[t] + citR[t]) atp[
                t] asp[t] ppi[t] amp[t] + 
            k3e2 k5e2 k8e2 k10e2 k12e2 atp[t] asp[t] ppi[t] amp[
                t] (as[t] + asR[t]));

Argininosuccinate lyase:

(*Velocity of the FORWARD argininosuccinate lyase reaction when the argininosuccinate is NOT radioactive.*)   vel3FOR := 
          aslAdj*(k1e3 k3e3 k5e3 as[t])*e30/denom3;

vel3REV := 
          aslAdj*(k2e3 k4e3 k6e3 fum[t] arg[t])*
            e30/denom3;

vel3FORRa := 
          aslAdj*(k1e3 k3e3 k5e3 asR[t])*e30/denom3;

vel3REVRa := 
          aslAdj*(k2e3 k4e3 k6e3 fum[t] argR[t])*
            e30/denom3;

denom3 := (k5e3 (k2e3 + k3e3) + (as[t] + 
                  asR[t]) k1e3 (k3e3 + k5e3) + 
            fum[t] k2e3 k4e3 + (arg[t] + argR[t] + 
                  fm*argQ[t]) k6e3 (k2e3 + k3e3) + 
            fum[t] (arg[t] + argR[t] + 
                  fm*argQ[t]) k4e3 k6e3 + (as[t] + asR[t]) fum[
                t] k1e3 k4e3);

Arginase:

argfac = argfac;

\(\(\( (*Rate\ equations\ for\ hydrolysis\ of\ arginine\ by\ arginase\ when\ the\ arginase\ is\ in\ the\ CYTOPLASM*) \)\(
\)\(
\)\(vel4FOR := 
      argfac*k1e4\ k3e4\ k4e4\ arg[t]\ e40/denom4; 
    vel4FORQ := argfac*fm\ k1e4\ k3e4\ k4e4\ argQ[t]\ e40/denom4;

    \(vel4FORRa := argfac*k1e4\ k3e4\ k4e4\ argR[t]\ e40/denom4;\)

    denom4 := 
      k4e4\ \((k2e4 + k3e4)\) + k5e4\ \((k2e4 + k3e4)\)\ orn[t] + 
        k1e4\ \((k3e4 + k4e4)\)\ \((arg[t] + argR[t] + 
              fm*argQ[t])\)\)\)\)

\(\(\( (*Equations\ for\ hydrolysis\ of\ arginine\ by\ arginase\ when\ the\ arginase\ is\ in\ the\ extracellular\ medium*) \)\(
\)\(
\)\( (*the\ denominator\ for\ extracellular\ arginine\ hydrolysis*) \)\(
\)\(denom4E := 
      k4e4\ \((k2e4 + k3e4)\) + k5e4\ \((k2e4 + k3e4)\)\ ornE[t] + 
        k1e4\ \((k3e4 + k4e4)\)\ \((argE[t] + argRE[t] + 
              fmE*argQE[
                  t])\)

 (*hydrolysis\ of\ extracellular\ arginine\ by\ extracellular\ arginase*) 

    vel4FORE := argfac*k1e4\ k3e4\ k4e4\ argE[t]\ e40E/denom4E; 
    vel4FORQE := argfac*fm\ k1e4\ k3e4\ k4e4\ argQE[t]\ e40E/denom4E;

    \(vel4FORRaE := argfac*k1e4\ k3e4\ k4e4\ argRE[t]\ e40E/denom4E;\)

    \)\)\)

Rate equations for appearance of various metabolites from their 'pools' or disappearance to their 'sinks'

vel5 := katp atppool

(*Pyrophosphate efflux to the XX sink.*) 
vel6 := kppi ppi[
            t]

vel7 := kfum fum[
            t]

vel8 := kcp cppool

vel108 := kcpei cpe[t]
vel109 := 
      kcpie cpims[t]

vel208 := kcpim cpims[t]
vel209 := 
      kcpmi cp[t]

(*RADIOACTIVE Carbamoyl phosphate exchange across the outer mitochondrial membrane.*)
\(vel108R := kcpei\ cpeR[t]

      vel109R := kcpie\ cpimsR[t]\)

(*RADIOACTIVE Carbamoyl phosphate exchange across the inner mitochondrial membrane.*)
\(vel208R := kcpim\ cpimsR[t]

      vel209R := kcpmi\ cpR[t]\)

vel9 := kasp asppool

vel10 := kamp amp[
            t]

vel11 := kpi pi[
            t]

vel12 := kure ure[t]

(*Radioactive urea efflux to the XX sink.*) 
vel12Ru := kure ureR[t]

Rate equations for membrane trasnport of ornithine and citrulline into and out of the mitochondria

vel14 := kornci orn[t]

vel15 := kornic ornims[t]

vel18 := kcitci cit[t]

      vel18Rc := kcitci citR[t]

vel19 := kcitic citims[t]

      vel19Rc := kcitic citimsR[t]

vel20 := kcitim citims[t]

      vel20Rc := kcitim citimsR[t]

vel21 := kcitmi citmat[t]

      vel21Rc := kcitmi citmatR[t]

vel22 := kornim ornims[t]

vel23 := kornmi ornmat[t]

vcitCE := 
      kcitCE cit[t]

    vcitRCE := kcitRCE citR[t]

    vatpCE := katpCE atp[t]

    vaspCE := kaspCE asp[t]

    vampCE := kampCE amp[t]

    vppiCE := kppiCE ppi[t]

    vasCE := kasCE as[t]

    vasRCE := kasRCE asR[t]

    vfumCE := kfumCE fum[t]

    vargQCE := kargQCE argQ[t]

    vargCE := kargCE arg[t]

    vargRCE := kargRCE argR[t]

    vureCE := kureCE ure[t]

    vureRCE := kureRCE ureR[t]

    vornCE := kornCE orn[t]

    vcpeCE := kcpeCE cpe[t]

    vcpeRCE := kcpeRCE cpeR[t]

(*the following are the rate equations for transport of metabolites from the extracellular medium to the cytoplasm.*) 

vcitCEREV := kcitCEREV citE[t]

    vcitRCEREV := kcitRCEREV citRE[t]

    vatpCEREV := katpCEREV atpE[t]

    vaspCEREV := kaspCEREV aspE[t]

    vampCEREV := kampCEREV ampE[t]

    vppiCEREV := kppiCEREV ppiE[t]

    vasCEREV := kasCEREV asE[t]

    vasRCEREV := kasRCEREV asRE[t]

    vfumCEREV := kfumCEREV fumE[t]

    vargQCEREV := kargQCEREV argQE[t]

    vargCEREV := kargCEREV argE[t]

    vargRCEREV := kargRCEREV argRE[t]

    vureCEREV := kureCEREV ureE[t]

    vureRCEREV := kureRCEREV ureRE[t]

    vornCEREV := kornCEREV ornE[t]

    vcpeCEREV := kcpeCEREV cpeE[t]

    vcpeRCEREV := kcpeRCEREV cpeRE[t]

Reactions

\(v[otc] := vel1FOR*\ matV - vel1REV*matV

      v[otcR] := vel1FORRc*matV - vel1REVRc*matV\)

\(v[citmi] := vel21*matV - vel20*imsV

      v[citmiR] := vel21Rc*matV - vel20Rc*imsV\)

\(v[citic] := \ vel19*imsV - vel18*cytV

      v[citicR] := vel19Rc*imsV - vel18Rc*cytV\)

\(v[ass] := vel2FOR*cytV - vel2REV*cytV

      v[assR] := vel2FORRc*cytV - vel2REVRa*cytV\)

\(v[asl] := \ vel3FOR*cytV - vel3REV*cytV

      v[aslR] := vel3FORRa*cytV - vel3REVRa*cytV\)

\(v[a] := \ vel4FOR*cytV

      v[aR] := vel4FORRa*cytV

      v[aQ] := \ vel4FORQ*cytV\)

\(v[aE] := \ vel4FORE*eV

      v[aRE] := vel4FORRaE*eV

      v[aQE] := \ vel4FORQE*eV\)

v[ornci] := 
        vel14*cytV - vel15*imsV

v[ornim] := 
        vel22*imsV - vel23*matV

v[cpin] := vel8*eV

v[cpci] := 
        vel108*cytV - vel109*imsV

v[cpim] := 
        vel208*imsV - vel209*matV

v[cpciR] := 
        vel108R*cytV - vel109R*imsV

v[cpimR] := 
        vel208R*imsV - vel209R*matV

v[atpin] := vel5*cytV

v[aspin] := vel9*cytV

v[ppisink] := vel6*cytV

v[fumsink] := vel7*cytV

v[pisink] := vel11*matV

v[ampsink] := vel10*cytV

v[usink] := vel12*cytV

      v[usinkR] := 0*vel12Ru*cytV

(*cytoplasmic membrane exchange of citrulline*) 
v[citCE] := 
      vcitCE*cytV - 
        vcitCEREV*
          eV
 (*cytoplasmic membrane exchange of radioactive citrulline*) 

    v[citRCE] := 
      vcitRCE*cytV - 
        vcitRCEREV*
          eV
 (*cytoplasmic membrane exchange of ATP*) 

    v[atpCE] := 
      vatpCE*cytV - 
        vatpCEREV*
          eV
 (*cytoplasmic membrane exchange of aspartate*) 

    v[aspCE] := 
      vaspCE*cytV - 
        vaspCEREV*
          eV
 (*cytoplasmic membrane exchange of AMP*) 

    v[ampCE] := 
      vampCE*cytV - 
        vampCEREV*
          eV
 (*cytoplasmic membrane exchange of pyrophosphate*) 

    v[ppiCE] := 
      vppiCE*cytV - 
        vppiCEREV*
          eV
 (*cytoplasmic membrane exchange of argininosuccinate*) 

    v[asCE] := 
      vasCE*cytV - 
        vasCEREV*
          eV
 (*cytoplasmic membrane exchange of radioactive argininosuccinate*) 

    v[asRCE] := 
      vasRCE*cytV - 
        vasRCEREV*
          eV
 (*cytoplasmic membrane exchange of fumarate*) 

    v[fumCE] := 
      vfumCE*cytV - 
        vfumCEREV*
          eV
 (*cytoplasmic membrane exchange of exogenous arginine*) 

    v[argQCE] := 
      vargQCE*cytV - 
        vargQCEREV*
          eV
 (*cytoplasmic membrane exchange of endogenous arginine*) 

    v[argCE] := 
      vargCE*cytV - 
        vargCEREV*
          eV
 (*cytoplasmic membrane exchange of radioactive arginine*) 

    v[argRCE] := 
      vargRCE*cytV - 
        vargRCEREV*
          eV
 (*cytoplasmic membrane exchange of urea*) 

    v[ureCE] := 
      vureCE*cytV - 
        vureCEREV*
          eV
 (*cytoplasmic membrane exchange of radioactive urea*) 

    v[ureRCE] := 
      vureRCE*cytV - 
        vureRCEREV*
          eV
 (*cytoplasmic membrane exchange of ornithine*) 

    v[ornCE] := 
      vornCE*cytV - 
        vornCEREV*
          eV
 (*cytoplasmic membrane exchange of carbamoylphosphate*) 

    v[cpeCE] := 
      vcpeCE*cytV - 
        vcpeCEREV*
          eV
 (*cytoplasmic membrane exchange of radioactive carbamoyl phosphate*) 

    v[cpeRCE] := vcpeRCE*cytV - vcpeRCEREV*eV

Equations, Matrices and NDSolveMatrix

The stiochiometry of each chemical and transport reaction is described here:

eqns = {{otc, 
            cp[t]/matV + ornmat[t]/matV -> 
              citmat[t]/matV + pi[t]/matV}, {otcR, 
            cpR[t]/matV + ornmat[t]/matV -> 
              citmatR[t]/matV + pi[t]/matV}, {citmi, 
            citmat[t]/matV -> citims[t]/imsV}, {citmiR, 
            citmatR[t]/matV -> citimsR[t]/imsV}, {citic, 
            citims[t]/imsV -> cit[t]/cytV}, {citicR, 
            citimsR[t]/imsV -> citR[t]/cytV}, {ass, 
            cit[t]/cytV + atp[t]/cytV + asp[t]/cytV -> 
              as[t]/cytV + amp[t]/cytV + ppi[t]/cytV}, {assR, 
            citR[t]/cytV + atp[t]/cytV + asp[t]/cytV -> 
              asR[t]/cytV + amp[t]/cytV + ppi[t]/cytV}, {asl, 
            as[t]/cytV -> fum[t]/cytV + arg[t]/cytV}, {aslR, 
            asR[t]/cytV -> fum[t]/cytV + argR[t]/cytV}, {a, 
            arg[t]/cytV -> orn[t]/cytV + ure[t]/cytV}, {aR, 
            argR[t]/cytV -> orn[t]/cytV + ureR[t]/cytV}, {aQ, 
            argQ[t]/cytV -> 
              orn[t]/cytV + ure[t]/cytV}, 
{aE, 
            argE[t]/eV -> ornE[t]/eV + ureE[t]/eV}, {aRE, 
            argRE[t]/eV -> ornE[t]/eV + ureRE[t]/eV}, {aQE, 
            argQE[t]/eV -> ornE[t]/eV + ureE[t]/eV}, {ornci, 
            orn[t]/cytV -> ornims[t]/imsV}, {ornim, 
            ornims[t]/imsV -> ornmat[t]/matV}, {cpin, 
            cppool/eV -> cpeE[t]/eV}, 
{cpci, 
            cpe[t]/cytV -> cpims[t]/imsV}, 
{cpim, 
            cpims[t]/imsV -> cp[t]/matV}, 
{cpciR, 
            cpeR[t]/cytV -> cpimsR[t]/imsV}, 
{cpimR, 
            cpimsR[t]/imsV -> cpR[t]/matV}, {atpin, 
            atppool/cytV -> atp[t]/cytV}, {aspin, 
            asppool/cytV -> asp[t]/cytV}, {ppisink, 
            ppi[t]/cytV -> XX}, {fumsink, 
            fum[t]/cytV -> XX}, {pisink, 
            pi[t]/matV -> XX}, {ampsink, 
            amp[t]/cytV -> XX}, {usink, 
            ure[t]/cytV -> XX}, {usinkR, 
            ureR[t]/cytV -> XX}, {citCE, 
            cit[t]/cytV -> citE[t]/eV}, {citRCE, 
            citR[t]/cytV -> citRE[t]/eV}, {atpCE, 
            atp[t]/cytV -> atpE[t]/eV}, {aspCE, 
            asp[t]/cytV -> aspE[t]/eV}, {ampCE, 
            amp[t]/cytV -> ampE[t]/eV}, {ppiCE, 
            ppi[t]/cytV -> ppiE[t]/eV}, {asCE, 
            as[t]/cytV -> asE[t]/eV}, {asRCE, 
            asR[t]/cytV -> asRE[t]/eV}, {fumCE, 
            fum[t]/cytV -> fumE[t]/eV}, {argQCE, 
            argQ[t]/cytV -> argQE[t]/eV}, {argCE, 
            arg[t]/cytV -> argE[t]/eV}, {argRCE, 
            argR[t]/cytV -> argRE[t]/eV}, {ureCE, 
            ure[t]/cytV -> ureE[t]/eV}, {ureRCE, 
            ureR[t]/cytV -> ureRE[t]/eV}, {ornCE, 
            orn[t]/cytV -> ornE[t]/eV}, {cpeCE, 
            cpe[t]/cytV -> cpeE[t]/eV}, {cpeRCE, 
            cpeR[t]/cytV -> cpeRE[t]/eV}};

Define the stoichiometry matrix (NM), the substrate matrix (SM) and the velocity matrix (VM):

NM := 
        NMatrix[eqns, {cppool, atppool, asppool, 
            XX}];  
NM // MatrixForm;

    
SM := 
        SMatrix[eqns, {cppool, atppool, asppool, 
            XX}];
SM // MatrixForm;

    
VM := 
        VMatrix[eqns, {cppool, atppool, asppool, 
            XX}];

Solve the timecourse:

Clear[argAdd]

    argAdd = 1*10^-3;

otcAdj = 
        1;

    assAdj = 2;

    aslAdj = 1;

    argfac = 1;

(*free mixing coefficient*) 
Clear[
      fm]

    fm = fm;

    fmE = 
        1;  (*free mixing coefficient outside the cells*)

\(\(\( (*maximumtime\ for\ the\ simulation, 
      s*) \)\(
\)\(\(maxT = 
        300;\)

 (*this\ command\ solves\ \((and\ times)\)\ the\ time\ course\ using\ NDSolve*) 

    First[
      Timing[sol = 
          Table[NDSolveMatrix[SM, NM, VM, initialconditions, {t, 0, maxT}, 
              WorkingPrecision -> 15, AccuracyGoal -> Automatic, 
              PrecisionGoal -> Automatic, 
              StartingStepSize -> Automatic, 
              MaxSteps -> 2000], {fm, 0.0, 1.0, 
              0.1}]
]]\)\)\)

RadioactivityDistribution::usage="RadioactivityDistribution[{(Subscript[time, 1]), (Subscript[time, 2]),...,(Subscript[time, n])}] generates a table of values for the percentage of each ornithine-dependent labelled metabolite. It also prints the number of nmol/ml of urea produced and the cpm in urea";
Begin["Private`"];

\(RadioactivityDistribution[
          timeT1_List] := 
Module[{d1b, ureR60b, argR60b, 
            asR60b, citR60b}, 

          d1b = \((ureRE[timeT1]\ eV + ureR[timeT1]*cytV + 
                argR[timeT1]\ cytV + argRE[timeT1]*eV + citR[timeT1]\ cytV + 
                citRE[timeT1]*eV + citmatR[timeT1]\ matV + 
                citimsR[timeT1]\ imsV + asR[timeT1]\ cytV + 
                asRE[timeT1]*eV)\); \  (*the\ denominator*) 

          ureR60b = 
            Flatten[Evaluate[\((ureR[timeT1]\ cytV + ureRE[timeT1]*eV)\)*100/
                      d1b /. sol]]; 

          argR60b = 
            Flatten[Evaluate[\((argR[timeT1]\ cytV + argRE[timeT1]*eV)\)*100/
                      d1b /. sol]]; 

          citR60b = 
            Flatten[Evaluate[\((citR[timeT1]\ cytV + citimsR[timeT1]\ imsV + 
                        citmatR[timeT1]\ matV + citRE[timeT1]\ eV)\)*100/d1b /. 
                  sol]]; 

          asR60b = 
            Flatten[Evaluate[\((asR[timeT1]\ cytV + asRE[timeT1]*eV)\)*100/
                      d1b /. sol]]; 

radb = {Prepend[ureR60b, "urea %"], 
              Prepend[argR60b, "arg %"], Prepend[citR60b, "cit %"], 
              Prepend[asR60b, "as %"]}; 

          TableForm[radb, TableDirections -> Column]]\)

End[];

Radioactivity::usage="Radioactivity[{(Subscript[time, 1]), (Subscript[time, 2]),...,(Subscript[time, n])}] generates a table of values for the total radioactivity of each ornithine-dependent labelled metabolite.";


\(Radioactivity[
        timeT2_List] := 
Module[{ureRcpm, argRcpm, 
          citRcpm}, 
ureRcpm = 
          Flatten[Evaluate[
              0.5*\((ureR[timeT2]*cytV + ureRE[timeT2]*eV)\)*10^9*530/
                    0.0045333 /. sol]]; 
argRcpm = 
          Flatten[Evaluate[
              0.5*\((argR[timeT2]*cytV + argRE[timeT2]*eV)\)*10^9*530/
                    0.0045333 /. sol]]; 
citRcpm = 
          Flatten[Evaluate[
              0.5*\((citR[timeT2]*cytV + citRE[timeT2]*eV + 
                      citimsR[timeT2]*imsV + citmatR[timeT2]*matV)\)*10^9*530/
                    0.0045333 /. 
                sol]]; 


title = {{metabolite, cpm}, Prepend[ureRcpm, "urea"], 
            Prepend[argRcpm, "arg"], Prepend[citRcpm, "cit"]}; 

        TableForm[title, TableDirections -> Column]]\)

RadPlot::usage="RadPlot[timeT1], plots the % radioactivity in each metabolite; urea=red, arginine=green, citrulline=blue and argininosuccinate=black";

Begin["Private`"];

\(RadPlot[timeT1_] := 
        Module[{d1b, ureR60c, argR60c, asR60c, citR60c}, 
          d1b = \((ureRE[t]\ eV + ureR[t]*cytV + argR[t]\ cytV + 
                argRE[t]*eV + citR[t]\ cytV + citRE[t]*eV + 
                citmatR[t]\ matV + citimsR[t]\ imsV + asR[t]\ cytV + 
                asRE[t]*eV)\); \  (*the\ denominator*) 

          ureR60c = \((ureR[t]\ cytV + ureRE[t]*eV)\)*100/d1b; 

          argR60c = \((argR[t]\ cytV + argRE[t]*eV)\)*100/d1b; 

          citR60c = \((citR[t]\ cytV + citimsR[t]\ imsV + citmatR[t]\ matV + 
                    citRE[t]\ eV)\)*100/d1b /. sol; 

          asR60c = \((asR[t]\ cytV + asRE[t]*eV)\)*100/
                d1b; 

Plot[
            Evaluate[{ureR60c, argR60c, citR60c, asR60c} /. sol], {t, 0, 
              timeT1}, PlotRange -> All, 
            PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], 
                RGBColor[0, 0, 1], RGBColor[0, 0, 0]}, 
            AxesLabel -> {"Time, s", "% radioactivity"}]];\)
End[];

Table Distributions

rd = RadioactivityDistribution[{60}]

radcpm = Radioactivity[{60}]

argQplotr = 
        Plot[Evaluate[1000*argQ[t] /. sol], {t, 0, 60}, 
          PlotRange -> All, 
          AxesLabel -> {"Time, s", "Added arginine, mM"}, 
          ImageSize -> 288*2];