Supplementary Notebook for
The Geometry of Multisite Phosphorylation
by Arjun Manrai & Jeremy Gunawardena

This should be read in conjunction with the paper. There are two changes of notation from the paper. Instead of E and F for kinase and phosphatase, respectively, we use K and F because E is a reserved symbol in Mathematica. Instead of a superscript to indicate the enzyme in a rate constant we use Roman letters for the kinase and Greek letters for the phosphatase. The correspondence between the paper's notation and the notebook's notation is summarised in the table below.

({{Paper, Notebook}, {E, K}, {F, F}, {(Subscript[((a^E)), u]), (Subscript[a, u])}, {(Subscript[((b^E)), u]), (Subscript[b, u])}, {(Subscript[((c^E)), u,v]), (Subscript[c, u,v])}, {(Subscript[((a^F)), u]), (Subscript[α, u])}, {(Subscript[((b^F)), u]), (Subscript[β, u])}, {(Subscript[((c^F)), u,v]), (Subscript[γ, u,v])}})

Definition of some useful lists, rules and functions

Lists of state variables

vars = {(Subscript[KS, "00"]),(Subscript[KS, "01"]),(Subscript[KS, "10"]),(Subscript[FS, "01"]),(Subscript[FS, "10"]),(Subscript[FS, "11"]),K,F, (Subscript[S, "00"]),(Subscript[S, "01"]),(Subscript[S, "10"]),(Subscript[S, "11"])};
linearvars = {(Subscript[S, "00"]), (Subscript[S, "01"]),(Subscript[S, "10"]),(Subscript[S, "11"])};
quadvars = { (Subscript[S, "00"]^2),(Subscript[S, "01"]^2),(Subscript[S, "10"]^2),(Subscript[S, "11"]^2),(Subscript[S, "00"]) (Subscript[S, "01"]),(Subscript[S, "00"]) (Subscript[S, "10"]),(Subscript[S, "00"]) (Subscript[S, "11"]),(Subscript[S, "01"]) (Subscript[S, "10"]),(Subscript[S, "01"]) (Subscript[S, "11"]) ,(Subscript[S, "10"]) (Subscript[S, "11"])};

Replacement rules

kinaseDistributive = (Subscript[c, "00","11"])-> 0;
phosphataseDistributive = (Subscript[γ, "11","00"])-> 0;
eliminateEnzymeSubstrates={(Subscript[KS, "00"])->(K (Subscript[S, "00"]))(Subscript[a, "00"])/((Subscript[b, "00"])+(Subscript[c, "00","01"])+(Subscript[c, "00","10"])+(Subscript[c, "00","11"])),
(Subscript[KS, "01"])-> (K (Subscript[S, "01"]))(Subscript[a, "01"])/((Subscript[b, "01"])+(Subscript[c, "01","11"])),
(Subscript[KS, "10"])-> (K (Subscript[S, "10"]))(Subscript[a, "10"])/((Subscript[b, "10"])+(Subscript[c, "10","11"])),
(Subscript[FS, "01"])-> (F (Subscript[S, "01"]))(Subscript[α, "01"])/((Subscript[β, "01"])+ (Subscript[γ, "01","00"])),
(Subscript[FS, "10"])-> (F (Subscript[S, "10"]))(Subscript[α, "10"])/((Subscript[β, "10"])+ (Subscript[γ, "10","00"])),
(Subscript[FS, "11"])-> (F (Subscript[S, "11"]))(Subscript[α, "11"])/((Subscript[β, "11"])+ (Subscript[γ, "11","01"])+(Subscript[γ, "11","10"])+(Subscript[γ, "11","00"]))};
replaceKandF = K -> t F;

Functions for inspecting rational expressions

degrees[x] returns the degree of the numerator of x, followed by the degree of denominator of x
cn[x, j] returns the coefficient of (t^j) in the numerator of x
cd[x, j] returns the coefficient of (t^j) in the denominator of x
These should only be used when x is expressed as a fraction, using the built-in function Together

quadcoeffs[x] returns a 10 x 2 matrix in which the first entry on row i is the i-th quadratic monomial in the quadvars list (above) and the second entry is the corresponding coefficient in x.

degrees[x_]:={Exponent[Numerator[x],t],Exponent[Denominator[x],t]};
cn[x_,j_]:=Simplify[Coefficient[Numerator[x],t,j]];
cd[x_,j_]:=Simplify[Coefficient[Denominator[x],t,j]];
cq[x_,i_]:={Part[quadvars,i],Coefficient[x,Part[quadvars,i]]};
quadcoeffs[x_]:=MatrixForm[Table[cq[x,i],{i,1,10}]];

Definition of the basis rate functions

Rate functions for multi-site phosphorylation with two sites, processive kinase and processive phosphatase

(Subscript[f, 1])= -(Subscript[a, "00"])K (Subscript[S, "00"])+(Subscript[b, "00"])(Subscript[KS, "00"])+(Subscript[γ, "01","00"])(Subscript[FS, "01"])+(Subscript[γ, "10","00"]) (Subscript[FS, "10"])+(Subscript[γ, "11","00"])(Subscript[FS, "11"]);
(Subscript[f, 2])= -(Subscript[a, "01"])K (Subscript[S, "01"])+(Subscript[b, "01"])(Subscript[KS, "01"])+(Subscript[c, "00","01"])(Subscript[KS, "00"]) - (Subscript[α, "01"])F (Subscript[S, "01"])+ (Subscript[β, "01"])(Subscript[FS, "01"])+(Subscript[γ, "11","01"])(Subscript[FS, "11"]);
(Subscript[f, 3])= -(Subscript[a, "10"])K (Subscript[S, "10"])+(Subscript[b, "10"])(Subscript[KS, "10"])+(Subscript[c, "00","10"])(Subscript[KS, "00"]) - (Subscript[α, "10"])F (Subscript[S, "10"])+(Subscript[β, "10"])(Subscript[FS, "10"])+(Subscript[γ, "11","10"])(Subscript[FS, "11"]);
(Subscript[f, 4])= -(Subscript[α, "11"])F (Subscript[S, "11"])+(Subscript[β, "11"])(Subscript[FS, "11"])+(Subscript[c, "01","11"])(Subscript[KS, "01"])+(Subscript[c, "10","11"])(Subscript[KS, "10"])+(Subscript[c, "00","11"])(Subscript[KS, "00"]);
(Subscript[f, 5])=  (Subscript[a, "00"])K (Subscript[S, "00"])-((Subscript[b, "00"])+(Subscript[c, "00","01"])+(Subscript[c, "00","10"])+(Subscript[c, "00","11"]))(Subscript[KS, "00"]);
(Subscript[f, 6])=  (Subscript[a, "01"])K (Subscript[S, "01"])-((Subscript[b, "01"])+(Subscript[c, "01","11"]))(Subscript[KS, "01"]);
(Subscript[f, 7])=  (Subscript[a, "10"])K (Subscript[S, "10"])-((Subscript[b, "10"])+(Subscript[c, "10","11"]))(Subscript[KS, "10"]);
(Subscript[f, 8])=  (Subscript[α, "01"])F (Subscript[S, "01"])-((Subscript[β, "01"])+ (Subscript[γ, "01","00"]))(Subscript[FS, "01"]);
(Subscript[f, 9])=  (Subscript[α, "10"])F (Subscript[S, "10"])-((Subscript[β, "10"])+ (Subscript[γ, "10","00"]))(Subscript[FS, "10"]);
(Subscript[f, 10])=  (Subscript[α, "11"])F (Subscript[S, "11"])-((Subscript[β, "11"])+ (Subscript[γ, "11","01"])+(Subscript[γ, "11","10"])+(Subscript[γ, "11","00"]))(Subscript[FS, "11"]);
(Subscript[f, 11])= -(Subscript[a, "00"])K (Subscript[S, "00"])+((Subscript[b, "00"])+(Subscript[c, "00","01"])+(Subscript[c, "00","10"])+(Subscript[c, "00","11"]))(Subscript[KS, "00"]) -(Subscript[a, "01"])K (Subscript[S, "01"])+((Subscript[b, "01"])+(Subscript[c, "01","11"]))(Subscript[KS, "01"])-(Subscript[a, "10"])K (Subscript[S, "10"])+((Subscript[b, "10"])+(Subscript[c, "10","11"]))(Subscript[KS, "10"]);
(Subscript[f, 12])= -(Subscript[α, "01"])F (Subscript[S, "01"])+((Subscript[β, "01"])+ (Subscript[γ, "01","00"]))(Subscript[FS, "01"])-(Subscript[α, "10"])F (Subscript[S, "10"])+((Subscript[β, "10"])+ (Subscript[γ, "10","00"]))(Subscript[FS, "10"])-(Subscript[α, "11"])F (Subscript[S, "11"])+((Subscript[β, "11"])+ (Subscript[γ, "11","01"])+(Subscript[γ, "11","10"])+(Subscript[γ, "11","00"]))(Subscript[FS, "11"]);

Check conservation of total kinase, total phosphatase and total substrate, respectively

Simplify[(Subscript[f, 11])+(Subscript[f, 5])+(Subscript[f, 6])+(Subscript[f, 7])]
Simplify[(Subscript[f, 12])+(Subscript[f, 8])+(Subscript[f, 9])+(Subscript[f, 10])]
Simplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])+(Subscript[f, 7])+(Subscript[f, 8 ])+(Subscript[f, 9 ])+(Subscript[f, 10 ])]

Kinase processive, phosphatase processive case I
Grobner basis

Calculate the Grobner basis

basis = Table[(Subscript[f, i]),{i,1,12}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];

The number of polynomials in the Grobner basis turns out to be the same as the number in the original basis

Length[grbasis]

The polynomials in the basis can be inspected by evaluating gr[i] for i from 1 to 12

gr[i_]:= Simplify[Part[grbasis,i]];

gr[12]

Calculate the phospho-form invariants

The first three polynomials in the Grobner basis involve only the state variables  (Subscript[S, 00]),(Subscript[S, 01]),(Subscript[S, 10]),(Subscript[S, 11]) and F. They form a basis for the 7th elimination ideal. However, each of these three polynomials is divisible by F. Dividing by F gives three invariants, (Subscript[p, 1]),(Subscript[p, 3]),(Subscript[p, 3]), for π(W).

(Subscript[p, 1])=gr[1]/F
(Subscript[p, 2])=gr[2]/F
(Subscript[p, 3])=gr[3]/F

The structure of these polynomials may be more easily understood using quadcoeffs

quadcoeffs[(Subscript[p, 1])]
quadcoeffs[(Subscript[p, 2])]
quadcoeffs[(Subscript[p, 3])]

Kinase distributive, phosphatase distributive (DD) case - I
Grobner basis

Since a Grobner basis was calculated above for the PP case, it might be thought that a Grobner basis for the DD case can be obtained by setting (Subscript[c, 00,11])= (Subscript[γ, 11,00])= 0. However, in general, specialisation of parameter values does not commute with forming Grobner bases. To avoid any difficulties arising from this, we calculate the Grobner basis for the DD case directly from the original rate functions for the dynamical system. We show below that, indeed, the resulting invariants are not the same polynomials as are obtained by setting (Subscript[c, 00,11])= (Subscript[γ, 11,00])= 0 in (Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3]).

Calculate the Grobner basis

basisDD = Table[((Subscript[f, i])/.kinaseDistributive)/.phosphataseDistributive,{i,1,12}];
grbasisDD=GroebnerBasis[basisDD,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];

Length[grbasisDD]

grDD[i_]:= Simplify[Part[grbasisDD,i]];

grDD[12]

Calculate the phospho-form invariants

(Subscript[((p^("DD"))), i_]):=grDD[i]/F;
(Subscript[((p^("DD"))), 1])
(Subscript[((p^("DD"))), 2])
(Subscript[((p^("DD"))), 3])

quadcoeffs[(Subscript[((p^("DD"))), 1])]
quadcoeffs[(Subscript[((p^("DD"))), 2])]
quadcoeffs[(Subscript[((p^("DD"))), 3])]

Specialisation of parameters does not commute with Grobner bases

FullSimplify[(((Subscript[p, 1])/.kinaseDistributive)/.phosphataseDistributive) - (Subscript[((p^("DD"))), 1])]
FullSimplify[(((Subscript[p, 2])/.kinaseDistributive)/.phosphataseDistributive) - (Subscript[((p^("DD"))), 2])]
FullSimplify[(((Subscript[p, 3])/.kinaseDistributive)/.phosphataseDistributive) - (Subscript[((p^("DD"))), 3])]

Calculate the planarity invariant

A linear combination of (Subscript[((p^("DD"))), 1]) and(Subscript[((p^("DD"))), 2]) simplifies to only involve 4 of the quadratic monomials

invariant =Collect[FullSimplify[((Subscript[b, "00"])+(Subscript[c, "00","01"])+(Subscript[c, "00","10"]))(Subscript[((p^("DD"))), 1])+((Subscript[c, "00","01"])+(Subscript[c, "00","10"]))((Subscript[β, "01"])+(Subscript[γ, "01","00"]))(Subscript[((p^("DD"))), 2])],{(Subscript[S, "10"]^2),(Subscript[S, "01"]^2)}]

quadcoeffs[invariant]

The (Subscript[μ, i]) are the coefficients of the linear form in the new variables (Subscript[y, i]).They will be used below. Note that they are all positive, in the sense explained in Section 2.4 of the paper.

(Subscript[μ, 1])=Coefficient[invariant,(Subscript[S, "01"]^2)]
 (Subscript[μ, 2])=Coefficient[invariant,(Subscript[S, "01"]) (Subscript[S, "10"])]
 (Subscript[μ, 3])=Coefficient[invariant,(Subscript[S, "10"]^2)]
 (Subscript[μ, 4])=-Coefficient[invariant,(Subscript[S, "00"]) (Subscript[S, "11"])]

Kinase processive, phosphatase processive case - II
Rational parametrisation

Calculate the rational functions (Subscript[r, i])

newfList = Table[((Subscript[f, i])/.eliminateEnzymeSubstrates)/.replaceKandF,{i,1,4}];
reducedmatrix = RowReduce[Last[Normal[CoefficientArrays[newfList,linearvars]]]];
(Subscript[r, i_]):=-Together[FullSimplify[reducedmatrix[[i,4]]]];
Table[degrees[(Subscript[r, i])],{i,1,3}]

Check the positivity property of the (Subscript[r, i])

The coefficients of the numerator and denominator of each (Subscript[r, i]) are all positive, in the sense of Section 2.4.

cn[(Subscript[r, 1]),0]
cn[(Subscript[r, 1]),1]
cn[(Subscript[r, 1]),2]
cd[(Subscript[r, 1]),0]
cd[(Subscript[r, 1]),1]
cd[(Subscript[r, 1]),2]
cd[(Subscript[r, 1]),3]

cn[(Subscript[r, 2]),0]
cn[(Subscript[r, 2]),1]
cd[(Subscript[r, 2]),0]
cd[(Subscript[r, 2]),1]
cd[(Subscript[r, 2]),2]

cn[(Subscript[r, 3]),0]
cn[(Subscript[r, 3]),1]
cd[(Subscript[r, 3]),0]
cd[(Subscript[r, 3]),1]
cd[(Subscript[r, 3]),2]

Calculate the (Subscript[y, i]) and determine the information for Figure 2

(Subscript[y, 1])= Together[((((Subscript[r, 2])))^2)/(Subscript[r, 1])];
(Subscript[y, 2])= Together[(Subscript[r, 2])(Subscript[r, 3])/(Subscript[r, 1])];
(Subscript[y, 3])= Together[((((Subscript[r, 3])))^2)/(Subscript[r, 1])];
Table[degrees[(Subscript[y, i])],{i,1,3}]

In view of their definition, the (Subscript[y, i]) inherit positivity from the (Subscript[r, i]). Hence, once their degrees are known, it is only necessary to determine their lowest monomials.

cn[(Subscript[y, 1]),0]
cn[(Subscript[y, 1]),1]
cd[(Subscript[y, 1]),0]

cn[(Subscript[y, 2]),0]
cn[(Subscript[y, 2]),1]
cd[(Subscript[y, 2]),0]

cn[(Subscript[y, 3]),0]
cn[(Subscript[y, 3]),1]
cd[(Subscript[y, 3]),0]

Inspect the rational parametrisation of the planar discrepancy for Figure 2

Δ = Together[(Subscript[μ, 1])(Subscript[y, 1]) + (Subscript[μ, 2])(Subscript[y, 2]) + (Subscript[μ, 3])(Subscript[y, 3])- (Subscript[μ, 4])];
degrees[Δ]

In view of its definition, the denominator of Δ inherits positivity from the (Subscript[r, i]), although it is still necessary to find its lowest monomial for Figure 2. The numerator of Δ may have general coefficients and the lowest and highest monomials need to be identified and their signs checked.

cn[Δ ,0]
cn[Δ ,4]
cd[Δ ,0]

Kinase distributive, phosphatase processive (DP) case
Rational parametrisation

The rational function (Subscript[r, i]) were obtained by Gaussian elimination over ℝ(a). It can be checked that, for the particular matrix that is used which arises from the rate equations, this process commutes with setting either or both of (Subscript[c, 00,11]) or (Subscript[γ, 11,00]) to 0. Hence the rational functions corresponding to the DP, PD and DD cases can be obtained from the (Subscript[r, i]) by setting (Subscript[c, 00,11])=0, (Subscript[γ, 11,00])=0 and (Subscript[c, 00,11])=(Subscript[γ, 11,00])=0, respectively. Note that they thereby inherit positivity from the (Subscript[r, i]).

Calculate the rational functions (Subscript[(r^DP), i])

(Subscript[((r^("DP"))), i_]):=Together[((Subscript[r, i])/.kinaseDistributive)];
Table[degrees[(Subscript[((r^("DP"))), i])],{i,1,3}]

Calculate the (Subscript[(y^DP), i]) and determine the information for Figure 2

(Subscript[((y^("DP"))), 1])= Together[((((Subscript[((r^("DP"))), 2])))^2)/(Subscript[((r^("DP"))), 1])];
(Subscript[((y^("DP"))), 2])= Together[(Subscript[((r^("DP"))), 2])(Subscript[((r^("DP"))), 3])/(Subscript[((r^("DP"))), 1])];
(Subscript[((y^("DP"))), 3])= Together[((((Subscript[((r^("DP"))), 3])))^2)/(Subscript[((r^("DP"))), 1])];
Table[degrees[(Subscript[((y^("DP"))), i])],{i,1,3}]

cn[(Subscript[((y^("DP"))), 1]),0]
cd[(Subscript[((y^("DP"))), 1]),0]

Inspect the rational parametrisation of the planar discrepancy for Figure 2

ΔDP = Together[(Subscript[μ, 1])(Subscript[((y^("DP"))), 1]) + (Subscript[μ, 2])(Subscript[((y^("DP"))), 2]) + (Subscript[μ, 3])(Subscript[((y^("DP"))), 3])- (Subscript[μ, 4])];
degrees[ΔDP]

cn[ΔDP,0]
cn[ΔDP,2]
cd[ΔDP,0]

Kinase processive, phosphatase distributive (PD) case
Rational parametrisation

Calculate the rational functions (Subscript[(r^PD), i])

(Subscript[((r^("PD"))), i_]):=Together[((Subscript[r, i])/.phosphataseDistributive)];
Table[degrees[(Subscript[((r^("PD"))), i])],{i,1,3}]

Calculate the (Subscript[(y^PD), i]) and determine the information for Figure 2

(Subscript[((y^("PD"))), 1])= Together[((((Subscript[((r^("PD"))), 2])))^2)/(Subscript[((r^("PD"))), 1])];
(Subscript[((y^("PD"))), 2])= Together[(Subscript[((r^("PD"))), 2])(Subscript[((r^("PD"))), 3])/(Subscript[((r^("PD"))), 1])];
(Subscript[((y^("PD"))), 3])= Together[((((Subscript[((r^("PD"))), 3])))^2)/(Subscript[((r^("PD"))), 1])];
Table[degrees[(Subscript[((y^("PD"))), i])],{i,1,3}]

cn[(Subscript[((y^("PD"))), 1]),0]
cn[(Subscript[((y^("PD"))), 1]),1]
cd[(Subscript[((y^("PD"))), 1]),0]

Inspect the rational parametrisation of the planar discrepancy for Figure 2

ΔPD = Together[(Subscript[μ, 1])(Subscript[((y^("PD"))), 1]) + (Subscript[μ, 2])(Subscript[((y^("PD"))), 2]) + (Subscript[μ, 3])(Subscript[((y^("PD"))), 3])- (Subscript[μ, 4])];
degrees[ΔPD]

cn[ΔPD,0]
cn[ΔPD,2]
cd[ΔPD,0]

Kinase distributive, phosphatase distributive (DD) case - II
Rational parametrisation

Calculate the rational functions (Subscript[(r^DD), i])

(Subscript[((r^("DD"))), i_]):=Together[(((Subscript[r, i])/.kinaseDistributive)/.phosphataseDistributive)];
Table[degrees[(Subscript[((r^("DD"))), i])],{i,1,3}]

Calculate the (Subscript[(y^DD), i])

(Subscript[((y^("DD"))), 1])= Together[((((Subscript[((r^("DD"))), 2])))^2)/(Subscript[((r^("DD"))), 1])];
(Subscript[((y^("DD"))), 2])= Together[(Subscript[((r^("DD"))), 2])(Subscript[((r^("DD"))), 3])/(Subscript[((r^("DD"))), 1])];
(Subscript[((y^("DD"))), 3])= Together[((((Subscript[((r^("DD"))), 3])))^2)/(Subscript[((r^("DD"))), 1])];

Check that the planar discrepancy is zero in this case

ΔDD = Simplify[(Subscript[μ, 1])(Subscript[((y^("DD"))), 1]) + (Subscript[μ, 2])(Subscript[((y^("DD"))), 2]) + (Subscript[μ, 3])(Subscript[((y^("DD"))), 3])- (Subscript[μ, 4])]

Plotting the ((Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3])) curves

Define replacement rules for parameter values

These values were drawn randomly from the uniform distribution on [0.00,5.00].

SetParms = {(Subscript[a, "00"]) -> 3.62, (Subscript[a, "01"]) -> 3.68,(Subscript[a, "10"]) -> 1.32,(Subscript[α, "11"]) -> 3.3,(Subscript[α, "10"]) -> 3.69,(Subscript[α, "01"]) -> 1.36,(Subscript[b, "00"]) -> 1.92,(Subscript[b, "01"]) ->4.66,(Subscript[b, "10"]) ->3.12,(Subscript[β, "11"]) -> 1.43,(Subscript[β, "10"]) ->3.18,(Subscript[β, "01"]) ->4.06,(Subscript[c, "00","01"]) -> 0.34, (Subscript[c, "00","10"]) -> 4.49,(Subscript[c, "01","11"]) ->3.07,(Subscript[c, "10","11"]) ->3.4,(Subscript[γ, "11","01"]) -> 2.88,(Subscript[γ, "11","10"]) -> 4.8,(Subscript[γ, "10","00"]) -> 0.85,(Subscript[γ, "01","00"]) -> 2.02};
SetDDParms = {(Subscript[c, "00","11"]) -> 0,(Subscript[γ, "11","00"]) -> 0};
SetDPParms = {(Subscript[c, "00","11"]) -> 0,(Subscript[γ, "11","00"]) -> 2.57};
SetPDParms = {(Subscript[c, "00","11"]) -> 4.83,(Subscript[γ, "11","00"]) -> 0};
SetPPParms = {(Subscript[c, "00","11"]) -> 4.83,(Subscript[γ, "11","00"]) -> 2.57};

Re-calculate the rational functions with these parameter values

(Subscript[ry, 1]):= ((Subscript[y, 1])/.SetParms)/.SetPPParms;
(Subscript[ry, 2]):= ((Subscript[y, 2])/.SetParms)/.SetPPParms;
(Subscript[ry, 3]):= ((Subscript[y, 3])/.SetParms)/.SetPPParms;
(Subscript[((ry^("DP"))), 1]):= ((Subscript[((y^("DP"))), 1])/.SetParms)/.SetDPParms;
(Subscript[((ry^("DP"))), 2]):= ((Subscript[((y^("DP"))), 2])/.SetParms)/.SetDPParms;
(Subscript[((ry^("DP"))), 3]):= ((Subscript[((y^("DP"))), 3])/.SetParms)/.SetDPParms;
(Subscript[((ry^("PD"))), 1]):= ((Subscript[((y^("PD"))), 1])/.SetParms)/.SetPDParms;
(Subscript[((ry^("PD"))), 2]):= ((Subscript[((y^("PD"))), 2])/.SetParms)/.SetPDParms;
(Subscript[((ry^("PD"))), 3]):= ((Subscript[((y^("PD"))), 3])/.SetParms)/.SetPDParms;
(Subscript[((ry^("DD"))), 1]):= ((Subscript[((y^("DD"))), 1])/.SetParms)/.SetDDParms;
(Subscript[((ry^("DD"))), 2]):= ((Subscript[((y^("DD"))), 2])/.SetParms)/.SetDDParms;
(Subscript[((ry^("DD"))), 3]):= ((Subscript[((y^("DD"))), 3])/.SetParms)/.SetDDParms;
plane = ((Subscript[μ, 1])x + (Subscript[μ, 2]) y + (Subscript[μ, 3])z- (Subscript[μ, 4])) /.SetParms;

The individual coordinate functions can be inspected by plotting them as a function of t

Table[Plot[(Subscript[ry, u]),{t,0,10}],{u,{1,2,3}}]
Table[Plot[(Subscript[((ry^("DP"))), u]),{t,0,2}],{u,{1,2,3}}]
Table[Plot[(Subscript[((ry^("PD"))), u]),{t,0,2}],{u,{1,2,3}}]
Table[Plot[(Subscript[((ry^("DD"))), u]),{t,0,2}],{u,{1,2,3}}]

Plot the (Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3]) curves in (ℝ^3) for all four cases

Show[ContourPlot3D[plane == 0 ,{x,0,0.8},{y,0,2.5},{z,0,11}, Mesh->None, ContourStyle -> Opacity[0.4], AxesLabel -> {"(Subscript[y, 1])","(Subscript[y, 2])","(Subscript[y, 3])"}, LabelStyle ->  Directive[FontFamily-> "Helvetica"]],ParametricPlot3D[{(Subscript[ry, 1]),(Subscript[ry, 2]),(Subscript[ry, 3])},{t,0.01,100},PlotStyle -> Directive[Purple,Thick],PlotRange -> Automatic],
ParametricPlot3D[{(Subscript[((ry^("DP"))), 1]),(Subscript[((ry^("DP"))), 2]),(Subscript[((ry^("DP"))), 3])},{t,0.01,100},PlotStyle -> Directive[Cyan,Thick],PlotRange -> Automatic],
ParametricPlot3D[{(Subscript[((ry^("PD"))), 1]),(Subscript[((ry^("PD"))), 2]),(Subscript[((ry^("PD"))), 3])},{t,0.01,100},PlotStyle -> Directive[Red,Thick],PlotRange -> Automatic],ParametricPlot3D[{(Subscript[((ry^("DD"))), 1]),(Subscript[((ry^("DD"))), 2]),(Subscript[((ry^("DD"))), 3])},{t,0.01,100},PlotStyle -> Directive[Blue,Thick],PlotRange -> Automatic]]