Add finite state projection example

[TorkelE]

An example of how to simulate the CME using FiniteStateProjection.jl.

I trialled having a collapsed "quick-start" example at the beginning (in the same way we had the environment set up in the CRN example library doc page).

[isaacsas]

Should we add a collapsible "installation instructions" too?

Have you double checked you have the axes correct on the heatmaps? (Sometimes some libraries flip the x and y axes I recall.)

It might be nice to also make the SS heatmap by running SSAs? That would be a nice comparison to show one gets the same result.

Finally, you probably want to pass vectors that label the X and Y ticks in the heatmaps, as otherwise I think they might start plotting at 1 and not zero. (Maybe plot the initial condition to confirm it looks ok and is labeled right?)

[isaacsas]

Suggested change

	Previously, we have shown how [*stochastic chemical kinetics*](@ref math_models_in_catalyst_sck_jumps) describe how chemical reaction network models can be [exactly simulated](@ref simulation_intro_jumps) (using e.g. [Gillespie's algorithm](https://en.wikipedia.org/wiki/Gillespie_algorithm)). We also described how the [SDE](@ref math_models_in_catalyst_cle_sdes) and [ODE](@ref math_models_in_catalyst_rre_odes) approaches were approximations of these jump simulations, and only valid for large copy numbers. To gain a good understanding of the system's time development, we typically have to carry out a large number of jump simulations. An alternative approach, however, is to instead simulate the *full probability distribution of the system*. This corresponds to the distribution from which these jump simulations are drawn.
	As previously discussed, [*stochastic chemical kinetics*](@ref math_models_in_catalyst_sck_jumps) models are mathematically given by jump processes that capture the exact times at which individual reactions occur, and the exact (integer) amounts of each chemical species at a given time. They represent a more microscopic model than [Chemical Langevin SDE](@ref math_models_in_catalyst_cle_sdes) models and [Reaction Rate Equation ODE](@ref math_models_in_catalyst_rre_odes) models, which can be interpreted as approximations to stochastic chemical kinetics models in the large population limit.
	One can study the dynamics of stochastic chemical kinetics models by simulating the stochastic processes using Monte Carlo methods. For example, they can be [exactly sampled](@ref simulation_intro_jumps) using [Stochastic Simulation Algorithms](https://en.wikipedia.org/wiki/Gillespie_algorithm) (SSAs), which are also often referred to as Gillespie's method. To gain a good understanding of a system's dynamics, one typically has to carry out a large number of jump process simulations to minimize sampling error. To avoid such sampling error, an alternative approach is to solve ODEs for the *full probability distribution* that these processes have a given value at each time. Knowing this distribution, one can then calculate any statistic of interest that can be sampled via running many SSA simulations.

[isaacsas]

Suggested change

[*The chemical master equation*](https://en.wikipedia.org/wiki/Master_equation) (CME) describes the time development of this distribution[^1]. In fact, this equation is at the core of chemical reaction network kinetics, with all other approaches (such as ODE, SDE, and Jump simulations) being derived as various approximations of it. The CME is a system of ODEs, with one variable *for each possible state of the system*. Each of the ODE's variables describes (the rate of change in) the probability of the system being in that state. For a system with a single species $X$, the CME looks like

[*The chemical master equation*](https://en.wikipedia.org/wiki/Master_equation) (CME) describes the time development of this probability distribution[^1], and is given by a (possibly infinite) coupled system of ODEs (with one ODE for each possible chemical state, i.e. number configuration, of the system). For a system with a single species $X$, the CME looks like

[isaacsas]

This is just wrong. In this notation the ODE rhs are all identical... I'd suggest just writing the general chemical master equation as in the math model description (or using that form but replacing vector X/x by scalars).

[isaacsas]

Suggested change

We note that, since (for almost all chemical reaction networks) there is a non-zero probability that $X$ is at any specific integer value, the CME have an infinite number of variables (corresponding to the infinite number of potential states). Hence, it cannot be solved practically. However, for high enough species values, the probability of the system attaining such values becomes negligibly small. Here, a truncated version of the CME can be solved practically. An approach for this is the *finite state projection*[^2]. Below we describe how to use the [FiniteStateProjection.jl](https://github.com/SciML/FiniteStateProjection.jl) package to solve the truncated CME (with the package's [documentation](https://docs.sciml.ai/FiniteStateProjection/dev/) providing a more extensive description). While the CEM approach can be very powerful, we note that even for systems with few species, the truncated CME typically have too many states to be feasibly solved.

For chemical reaction networks in which the total population is bounded, the CME corresponds to a finite set of ODEs. In contrast, for networks in which the system can (in theory) become unbounded, such as networks that include zero order reactions like $\varnothing \to A$, the CME will correspond to an infinite set of ODEs. Even in the finite case, the number of ODEs corresponds to the number of possible state vectors, i.e. vectors with components representing the integer populations of each species in the network. Therefore, for even simple reaction networks there can be many more ODEs than can be represented in typical levels of computer memory, and it becomes infeasible to numerically solve the full system of ODEs that correspond to the CME. However, in many cases the probability of the system attaining species values outside some small range can become negligibly small. Here, a truncated, approximating, version of the CME can be solved practically. An approach for this is the *finite state projection method*[^2]. Below we describe how to use the [FiniteStateProjection.jl](https://github.com/SciML/FiniteStateProjection.jl) package to solve the truncated CME (with the package's [documentation](https://docs.sciml.ai/FiniteStateProjection/dev/) providing a more extensive description). While the CME approach can be very powerful, we note that even for systems with a few species, the truncated CME typically has too many states for it to be feasible to solve the full set of ODEs.

[isaacsas]

Suggested change

	As a first step, we import the FiniteStateProjection package. Next, we convert our [`ReactionSystem`](@ref) to a `FSPSystem` (from which we later will generate an ODE).
	As a first step, we import the FiniteStateProjection package. Next, we convert our [`ReactionSystem`](@ref) to a `FSPSystem` (from which we later will generate the ODEs that correspond to the truncated CME).