Supplemental Code for the Steady-State Reduction (SSR) and SSR-generated Parameter Change (SPARC) Methods
============================

This repository stores supplemental code for the papers:

+ "Steady State Reduction of generalized Lotka-Volterra systems in the
  microbiome," by Eric Jones and Jean Carlson, published in Physical Review E
  in 2019
  [(link)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.99.032403),
  (called the SSR paper) and
+ "Control of ecological outcomes through deliberate parameter changes in a
  model of the gut microbiome", by Zipeng Wang, Eric Jones, Joshua Mueller, and
  Jean Carlson, accepted to Physical Review E in 2020 [(arXiv
  link)](https://arxiv.org/abs/1912.03412) (called the SPARC paper).

The code in `steady_state_reduction_example.py` generates Fig. 2 of the SSR
paper.  Steady State Reduction (SSR) approximates a high-dimensional gLV system
by two-dimensional (2D) subspaces that each obey gLV dynamics and span a pair
of high-dimensional steady states. Therefore, the process of SSR determines the
2D parameters that best approximate in-plane dynamics of the original
high-dimensional system.

The code in `SPARC_example.py` generates Fig. 3 of the SPARC paper. The control
framework SPARC (SSR-generated Parameter Change) alters the steady state
outcome in bistable regions of high-dimensional gLV systems. SPARC operates by
deliberately changing interaction parameters of the high-dimensional model:
once the parameter change is made, the dynamics of an initial condition that
was initially flowing towards some steady state A will be altered so that it
now flows towards some steady state B.  Here it is applied to an 11-dimensional
generalized Lotka-Volterra (gLV) model fit to gut microbial data from a C.
difficile infection (CDI) mouse experiment.

`steady_state_reduction_example.py`
-----------------------------------

The structure of this code is, briefly:

1) Import the parameters that characterize a high-dimensional system (growth
   rates \rho and interactions K), and store them in the class `Params`. In this
   code, as in the paper, we use the 11-dimensional parameters fit by Stein et al.
   (PLOS Comp Bio, 2013). To apply SSR to your own high-dimensional system, create
   your own `Params` class specific to your data.

2) Generate a pair of steady states to use in SSR. This code finds steady
   states of the high-dimensional model with `get_all_stein_steady_states`, but to
   apply SSR to a system you are interested in, identify the pair of
   high-dimensional steady states of your particular system you would like to
   study.
3) Generate the 2D parameters with SSR. Parameters of this 2D system are stored
   in the class `ssrParams`; generate this class using the parameters and steady
   state of your interest.
4) Plot and compare (in-plane) trajectories in the high-dimensional and
   SSR-reduced systems, saved in `SSR_demo_1.pdf`.
5) Plot the separatrix of the SSR-reduced 2D system, saved in `SSR_demo_2.pdf`.
6) Plot the in-plane separatrix of the high-dimensional system, saved in
   `SSR_demo_2.pdf`. Depending on the spatial resolution of the generated
   high-dimensional in-plane separatrix, this simulation can take a while-- use
   the `load_data` parameter to save a generated separatrix to the disk.

To use SSR yourself, generate your own `Params` class, say `p`, and generate the pair of
steady states you intend to use for SSR, say `ssa` and `ssb`. Then generate new
ssrParams with `s = ssrParams(p, ssa, ssb)`. You can check the validity of SSR
as implemented in your own system with the three `plot_*` commands.

`SPARC_example.py`
-----------------------------------

In total, this code:

1) imports the SSR formulae by Jones and Carlson (Phys. Rev. E 2019)
   (referred to as the SSR paper) and gLV system parameters fit by Stein et al.
   (PLOS Comp. Bio. 2013);
2) compresses the high-dimensional parameters into 2D SSR-generated
   parameters and generate trajectories in both systems starting at the initial
   condition (0.5,0.5);
3) finds a 2D parameter change that alters the steady-state outcome of this
   initial condition and plots the resulting new trajectory in the reduced 2D
   system; and
4) generates a high-dimensional parameter modification that corresponds to
   the 2D parameter change and plots the corresponding trajectory.

To apply the SPARC method to your own gLV model, in the main function at the
end of the file, create your own `Params` container class with your own growth
rates `rho`, interaction matrix `K`, and labels `labels` with the command `my_p
= ssr.Params([labels, rho, K])`. Name the two steady states that define your
bistable region of interest `ssa` and `ssb`, and compute the SSR-reduced
parameter set with the command `my_s = ssr.ssrParams(my_p, ssa, ssb)`. Finally, to
generate the four-panel plot for your own gLV model starting from an initial
condition `my_IC`, use the command `plot_all_panels(my_p, my_s, my_IC)`.

Please email ewj@physics.ucsb.edu with questions or bugs.
