EXAMPLE 1: Convert a probabilistic population protocol to an ODE system using the balancing equation: (See Huang & Huls, 2022 ) and (Bournez et. al, 2012)¶
If $P$ is a system of reactions of the form $q_i q_j \rightarrow \alpha_{ijkl} q_k q_l$, and for each pair $(qi qj)$ the sum of all $\alpha_{ij \cdot \cdot} = 1$, then the equations that govern the entire system taken as an ODE are $$b(\bf{x}) = \sum_{i,j} ((x_i * x_j) (-e_i -e_j + \sum_{k,l} \alpha_{i,j,k,l}(e_k + e_l))).$$
Here,
- $e_n$ is the vector with a 1 in the $n$ index, and $0$ for all others.
- The outer sum is taken over all the different pairs of reactants.
- The inner sum is taken over all the different pairs of products, given some pair of reactants.
Suppose we have the following system $P$:
$X + Y \overset{1/3}{\rightarrow} X + C$
$X + Y \overset{2/3}{\rightarrow} X + D$
$C + D \overset{1/4}{\rightarrow} C + D$
$C + D \overset{3/4}{\rightarrow} C + Y$
Note that:
There is a certain subtle property here: For all reactions that have A and B as reactants, if one of these is a catalyst in some such reaction, then the other is not a catalyst in any of those reactions. For example, to maintain this property we could not include (even if we had probability left over) a reaction like $X + Y \overset{...}{\rightarrow} D + Y$ along with the above reactions. This is called a one-way protocol.
i.e. "Every interaction changes at most the responder's state; thus it can be implemented with one-way communication." (Angluin, Aspnes, Eisenstat - A Simple Population Protocol for Fast Robust Approximate Majority)
The system is symmetric in that (X + Y) in the reactants means the same thing as (Y+X) in the reactants. That is, it has a symmetric interpretation. We could just as well interpret it to be asymmetric and note there just happen to be no (Y+X) reactions.
This system is probabilistic. When X and Y run into each other, there is a 1/3 chance that X and C are produced, i.e. that $Y$ is replaced by $C$.
$\alpha_{x,y,x,c} = 1/3$ and $\sum \alpha_{x,y,\cdot,\cdot} = 1$.
If we solve for $b$, we get
$b(\bf{x}) = $
$xy * \begin{bmatrix} -1 \\ -1 \\ 0 \\ 0 \end{bmatrix} + xy \begin{bmatrix} 1/3 \\ 0 \\ 1/3 \\ 0 \end{bmatrix} + xy \begin{bmatrix} 2/3 \\ 0 \\ 0 \\ 2/3 \end{bmatrix}$
$+ cd \begin{bmatrix} 0 \\ 0 \\ -1 \\ -1 \end{bmatrix} + cd \begin{bmatrix} 0 \\ 0 \\ 1/4 \\ 1/4 \end{bmatrix} + cd \begin{bmatrix} 0 \\ 3/4 \\ 3/4 \\ 0 \end{bmatrix},$
which gives us the ODE system
$b(\bf{x}) = \begin{bmatrix} x' \\ y' \\ c' \\ d' \end{bmatrix} = \begin{bmatrix} 0 \\ 3/4 cd - xy \\ 1/3 xy \\ 2/3 xy -3/4 cd \end{bmatrix}$
Finally, note that we can solve for these values just like we would solve for ODEs of chemical networks if we choose to treat the probabilities (over each arrow) as rate constants. Check that the solution is indeed the same.