> #OK to post #Timothy Nasralla, HW21, 11/15/21
> read "/Users/tan88/OneDrive - Rutgers University/DMB.txt"
                          First Written: Nov. 2021 


This is DMB.txt, A Maple package to explore Dynamical models in Biology (both 

   discrete and continuous)


   accompanying the class Dynamical Models in Biology, Rutgers University. 

      Taught by Dr. Z. (Doron Zeilbeger) 



               The most current version is available on WWW at:

          http://sites.math.rutgers.edu/~zeilberg/tokhniot/DMB.txt .

           Please report all bugs to: DoronZeil at gmail dot com .


             For general help, and a list of the MAIN functions,

        type "Help();". For specific help type "Help(procedure_name);" 


                        ------------------------------

            For a list of the supporting functions type: Help1();

             For help with any of them type: Help(ProcedureName);


                        ------------------------------

  For a list of the functions that give examples of Discrete-time dynamical 

     systems (some famous), type: HelpDDM();


             For help with any of them type: Help(ProcedureName);


                        ------------------------------

 For a list of the functions continuous-time dynamical systems (some famous) 

    type: HelpCDM();


             For help with any of them type: Help(ProcedureName);


                        ------------------------------

> Help(ChemoStat)
 ChemoStat(N,C,a1,a2): The Chemostat continuous-time dynamical system with 

    N=Bacterial poplulation densitty, and C=nutient Concentration in growth 

    chamber (see Table 4.1 of Edelstein-Keshet, p. 122)


  with paramerts a1, a2, Equations (19a_, (19b) in Edelestein-Keshet p. 127 

     (section 4.5, where they are called alpha1, alpha2). a1 and a2 can be 

     symbolic or numeric. Try:


                            ChemoStat(N,C,a1,a2);

                             ChemoStat(N,C,2,3);

> Chemo1 := ChemoStat(N,C,2,3)

> TimeSeries(Chemo1, [N,C], [3, 3], 0.01, 10, 1); TimeSeries(Chemo1, [N,C], [3,
> 3], 0.01, 10, 2)


> PhaseDiag(Chemo1, [N,C], [3,3], 0.01, 10)

> SEquP(Chemo1, [N,C])
                                  {[4., 1.]}

> #As shown by the horizontal asymptotes approached by the respective graphs of
> N and C to time, the equilibrium point is [4,1], which was confirmed by the
> SEquP function.
> Chemo2 := ChemoStat(N,C,2,2)

> TimeSeries(Chemo2, [N,C], [2, 2], 0.01, 10, 1); TimeSeries(Chemo2, [N,C], [2,
> 2], 0.01, 10, 2)


> PhaseDiag(Chemo2, [N,C], [2,2], 0.01, 10)

> SEquP(Chemo2, [N,C])
                                  {[2., 1.]}

> #Similar to problem 1, the equilibrium point found that was stable is [2,1],
> as shown in the graphs and confirmed by SEquP.
> Chemo3 := ChemoStat(N,C,4,1)

> TimeSeries(Chemo3, [N,C], [1, 1], 0.01, 10, 1); TimeSeries(Chemo3, [N,C], [1,
> 1], 0.01, 10, 2)


> PhaseDiag(Chemo2, [N,C], [1,1], 0.01, 10)

> SEquP(Chemo3, [N,C])
                        {[2.666666667, 0.3333333333]}

> 

> Help(GeneNet)
> 
  GeneNet(a0,a,b,n,m1,m2,m3,p1,p2,p3): The contiuous-time dynamical system, 

     with quantities m1,m2,m3,p1,p2,p3, due to M. Elowitz and S. Leibler


  described in the Ellner-Guckenheimer book, Eq. (4.1) (chapter 4,  p. 112)

  and parameers a0 (called alpha_0 there),a (called alpha there), b (called 

     beta there) and n. Try:


                   GeneNet(0,0.5,0.2,2,m1,m2,m3,p1,p2,p3);

> #Note that since this transformation does not take place in R2, the PhaseDiag
> function wont work on it. Wow1 := GeneNet(1,2,3,1,m1,m2,m3,p1,p2,p3)

> TimeSeries(Wow1, [m1,m2,m3,p1,p2,p3], [1, 1,1,1,1, 1], 0.01, 10, 1);
> TimeSeries(Wow1, [m1,m2,m3,p1,p2,p3], [1, 1,1,1,1, 1], 0.01, 10, 2);
> TimeSeries(Wow1, [m1,m2,m3,p1,p2,p3], [1, 1,1,1,1, 1], 0.01, 10, 3);
> TimeSeries(Wow1, [m1,m2,m3,p1,p2,p3], [1, 1,1,1,1, 1], 0.01, 10, 4);
> TimeSeries(Wow1, [m1,m2,m3,p1,p2,p3], [1, 1,1,1,1, 1], 0.01, 10, 5);
> TimeSeries(Wow1, [m1,m2,m3,p1,p2,p3], [1, 1,1,1,1, 1], 0.01, 10, 6);






> SEquP(Wow1, [m1,m2,m3,p1,p2,p3])
{[1.732050808, 1.732050808, 1.732050808, 1.732050808, 1.732050808, 1.732050808]}

> #As shown by the phase planes that all approach the horizontal asymptote of
> 1.732 ish (sqrt 3), it is an equilibrium point.
> Wow2 := GeneNet(1,1,1,1,m1,m2,m3,p1,p2,p3)

> TimeSeries(Wow2, [m1,m2,m3,p1,p2,p3], [1, 2,1,2,1, 2], 0.01, 10, 1);
> TimeSeries(Wow2, [m1,m2,m3,p1,p2,p3], [1, 2,1,2,1, 2], 0.01, 10, 2);
> TimeSeries(Wow2, [m1,m2,m3,p1,p2,p3], [1, 2,1,2,1, 2], 0.01, 10, 3);
> TimeSeries(Wow2, [m1,m2,m3,p1,p2,p3], [1, 2,1,2,1, 2], 0.01, 10, 4);
> TimeSeries(Wow2, [m1,m2,m3,p1,p2,p3], [1, 2,1,2,1, 2], 0.01, 10, 5);
> TimeSeries(Wow2, [m1,m2,m3,p1,p2,p3], [1, 2,1,2,1, 2], 0.01, 10, 6);






> SEquP(Wow2, [m1,m2,m3,p1,p2,p3])
{[1.414213562, 1.414213562, 1.414213562, 1.414213562, 1.414213562, 1.414213562]}

> #In this situation, the GeneNet system approached the stable equilibrium of
> sqrt 2 for all values.
> Wow3 := GeneNet(2,4,6,0,m1,m2,m3,p1,p2,p3)

> TimeSeries(Wow3, [m1,m2,m3,p1,p2,p3], [1, 2,3,1,2, 3], 0.01, 10, 1);
> TimeSeries(Wow3, [m1,m2,m3,p1,p2,p3], [1, 2,3,1,2, 3], 0.01, 10, 2);
> TimeSeries(Wow3, [m1,m2,m3,p1,p2,p3], [1, 2,3,1,2, 3], 0.01, 10, 3);
> TimeSeries(Wow3, [m1,m2,m3,p1,p2,p3], [1, 2,3,1,2, 3], 0.01, 10, 4);
> TimeSeries(Wow3, [m1,m2,m3,p1,p2,p3], [1, 2,3,1,2, 3], 0.01, 10, 5);
> TimeSeries(Wow3, [m1,m2,m3,p1,p2,p3], [1, 2,3,1,2, 3], 0.01, 10, 6);






> SEquP(Wow3,[m1,m2,m3,p1,p2,p3])
                          {[4., 4., 4., 4., 4., 4.]}

> #All points moved towards 4, making it a stable equilibrium point.

> Help(Lotka)
Lotka(r1,k1,r2,k2,b12,b21,N1,N2): The Lotka-Volterra continuous-time 

   dynamical system, Eqs. (9a),(9b)  (p. 224, section 6.3) of Edelstein-Keshet


with popoluations N1, N2, and parameters r1,r2,k1,k2, b12, b21 (called there 

   beta_12 and beta_21)


                                     Try:

                      Lotka(r1,k1,r2,k2,b12,b21,N1,N2);

                          Lotka(1,2,2,3,1,2,N1,N2);

> #Since this function is rooted in R2, the PhaseDiag function will work on it.
> Lokta1 := Lotka(1,2,3,4,5,6,N1,N2)

> TimeSeries(Lokta1, [N1,N2], [1, 1], 0.01, 10, 1); TimeSeries(Lokta1, [N1,N2],
> [1, 1], 0.01, 10, 2);


> PhaseDiag(Lokta1, [N1,N2], [1, 1], 0.01, 10)

> SEquP(Lokta1, [N1,N2])
                             {[0., 4.], [2., 0.]}

> #The phase diagram does seemingly approach two values, [2,0] or [0,4]. That
> being said, the TimeSeries reflects the first stable equilibrium point, [0,4].
> Lokta2 := Lotka(1,1,1,1,1,1,N1,N2)

> TimeSeries(Lokta2, [N1,N2], [2, 2], 0.01, 10, 1); TimeSeries(Lokta2, [N1,N2],
> [2, 2], 0.01, 10, 2);


> PhaseDiag(Lokta2, [N1,N2], [2, 2], 0.01, 10)

> SEquP(Lokta2, [N1,N2])
                                      {}

> #Im somewhat surprised that the SEquP returned the empty set, since both time
> series and the PhaseDiagram reflect that the function approached (0.5,0.5).
> Lokta3 := Lotka(2,1,2,1,2,1,N1,N2)

> TimeSeries(Lokta3, [N1,N2], [1, 2], 0.01, 10, 1); TimeSeries(Lokta3, [N1,N2],
> [1, 2], 0.01, 10, 2);


> PhaseDiag(Lokta3, [N1,N2], [1, 2], 0.01, 10)

> SEquP(Lokta3, [N1,N2])
                                  {[0., 1.]}

> #As shown, this transformation approaches (0,1) as a stable equilibrium point.


> Help(Volterra)
    Volterra(a,b,c,d,x,y): The (simple, original) Volterra predator-prey 

       continuous-time dynamical system with parameters a,b,c,d


       Given by Eqs. (7a) (7b) in Edelstein-Keshet p. 219 (section 6.2).

                      a,b,c,d may be symbolic or numeric

                                    Try: 

                            Volterra(a,b,c,d,x,y);

                            Volterra(1,2,3,4,x,y);

> Volt1 := Volterra(1,2,3,4,x,y)

> TimeSeries(Volt1, [x,y], [1,1], 0.01, 10, 1); TimeSeries(Volt1, [x,y], [1,1],
> 0.01, 10, 2)


> PhaseDiag(Volt1, [x,y], [1,1], 0.01, 10)

> SEquP(Volt1, [x,y])
                                      {}

> #As seen in the phase diagram and the time series, both x and y oscillate
> which leads to no stable equilibrium points.

> Volt2 := Volterra(1,2,2,1,x,y)

> TimeSeries(Volt2, [x,y], [2,1], 0.01, 10, 1); TimeSeries(Volt2, [x,y], [2,1],
> 0.01, 10, 2)


> PhaseDiag(Volt2, [x,y], [2,1], 0.01, 10)

> SEquP(Volt2, [x,y])
                                      {}

> #Similar to the first problem, this system for Volt has no stable equilibrium
> points.
> Volt3 := Volterra(1,1,1,1,x,y)

> TimeSeries(Volt3, [x,y], [1,1], 0.01, 10, 1); TimeSeries(Volt3, [x,y], [1,1],
> 0.01, 10, 2)


> PhaseDiag(Volt3, [x,y], [1,1], 0.01, 10)

> SEquP(Volt3, [x,y])
                                      {}

> #Weirdly enough, the phase diagram reflects that [1,1] is a stable fixed point
> since it does not tend to anywhere, and both time series show it to be a
> horizontal asymptote; however the SEquP function returns no stable fixed
> points.

> Help(VolterraM)
VolterraM(a,b,c,d,x,K,y): The MODIFIED Volterra predator-prey continuous-time 

   dynamical system with parameters a,b,c,d,K


       Given by Eqs. (8a) (8b) in Edelstein-Keshet p. 220 (section 6.2).

                     a,b,c,d ,Kmay be symbolic or numeric

                                    Try: 

                          VolterraM(a,b,c,d,K,x,y);

                          VolterraM(1,2,3,4,3,x,y);

> VoltM1 := VolterraM(1,2,3,4,3,x,y)

> TimeSeries(VoltM1,[x,y],[1,.5],0.01,10,1);
> TimeSeries(VoltM1,[x,y],[1,.5],0.01,10,2)


> PhaseDiag(VoltM1,[x,y],[1,.5],0.01,10)

> SEquP(VoltM1,[x,y])
                             {[1., 0.3750000000]}

> #While the time series are oscillating, their oscillations seem to be
> decreasing implying an end, and the phase diagram is seemingly approaching the
> stable equilibrium found by the function SEquP
> 

> TimeSeries(VoltM2,[x,y],[0.7,0.3],0.01,10,1);
> TimeSeries(VoltM2,[x,y],[0.7,0.3],0.01,10,2);


> PhaseDiag(VoltM2,[x,y],[0.7,0.3],0.01,10)

> SEquP(VoltM2,[x,y])
                             {[1., 0.5000000000]}

> #As seen in the phase diagram and time series, the stable equilibrium point is
> (1,0.5)
> VoltM3 := VolterraM(0.1,0.2,0.3,0.4,0.5,x,y)

> TimeSeries(VoltM3,[x,y],[.1,10],0.01,10,1);
> TimeSeries(VoltM3,[x,y],[.1,10],0.01,10,2);


> PhaseDiag(VoltM3, [x,y], [0.1,10], 0.01,10)

> SEquP(VoltM3,[x,y])
                             {[0.4000000000, 0.]}

> #While the time series for x doesn't fully reach 0.4, it approaches it rapidly
> and the phase diagram reflects a very miniscule equilibrium point which was
> then specified by the SEquP function to be [0.4,0].