# The moments of the distribution for mRNA in a two- and three-state promoter.

(c) 2017 Manuel Razo. This work is licensed under a [Creative Commons Attribution License CC-BY 4.0](https://creativecommons.org/licenses/by/4.0/). All code contained herein is licensed under an [MIT license](https://opensource.org/licenses/MIT)


```python
import pickle
import dill # Library that we will use to export lambdify functions
# Our numerical workhorses
import numpy as np
# To compute symbolic expressions
import sympy
sympy.init_printing(use_unicode=True, use_latex=True) # print outputs in LaTeX

# Import matplotlib stuff for plotting
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib as mpl

# Seaborn, useful for graphics
import seaborn as sns

# Import the utils for this project
import chann_cap_utils as chann_cap

# Set PBoC plotting style
chann_cap.set_plotting_style()

# Magic function to make matplotlib inline; other style specs must come AFTER
%matplotlib inline

# This enables SVG graphics inline (only use with static plots (non-Bokeh))
%config InlineBackend.figure_format = 'svg'

# figdir = '../../fig/chemical_master_analytic_mRNA/'
tmpdir = '../../tmp/'
```

### $\LaTeX$ macros
$\newcommand{kpon}{k^p_{\text{on}}}$
$\newcommand{kpoff}{k^p_{\text{off}}}$
$\newcommand{kron}{k^r_{\text{on}}}$
$\newcommand{kroff}{k^r_{\text{off}}}$
$\newcommand{rm}{r _m}$
$\newcommand{gm}{\gamma _m}$
$\newcommand{mm}{\left\langle m \right\rangle}$

# Two-state promoter (mRNA)

In this section we will perform the linear algebra to solve for the steady-state moments of the mRNA distribution.

Up to this stage we wrote the chemical master equation in matrix notation and did some tricks when it came to compute the $n^{\text{th}}$ moment. Now we will use the equations we obtained for the time derivative of the mRNA distribution moments and set that to steady-state to solve for the moments themselves.

Let's first define the `sympy` variables that we will need for the two- and three-state promoter.


```python
# Define rate constant variables
kp_off, kp_on = sympy.symbols('{k_{off}^{(p)}} {k_{on}^{(p)}}')
kr_off, kr_on = sympy.symbols('{k_{off}^{(r)}} {k_{on}^{(r)}}')

# Define degradation rate and production rate
rm, gm = sympy.symbols('r_m gamma_m')

kp_off, kp_on, kr_off, kr_on, rm, gm
```

Let's now define a function that generates a `sympy` variable of the form $\left\langle x^n \right\rangle$ to easily define the variables to save the moments.


```python
def moment_sympy(var, moment):
    '''
    Shortcut to generate variables of the form <var**n> in sympy to facilitate the
    definition of the moment variables
    Parameters
    ----------
    var : string.
        variable name string that should be put inside the <> parenthesis.
    moment : int.
        exponent of the moment to be defined.
    '''
    if moment == None:
        return sympy.Symbol('{\\langle ' + var + '\\rangle}')
    else:
        return sympy.Symbol('{\\langle ' + var + '^{:d}'.format(moment) +\
                            '\\rangle}')
```

Now let's define the moment variables for the unregulated promoter.


```python
# Define the string for the states
m_unreg = ['m_{P_{unreg}}', 'm_{E_{unreg}}']

# Define the zeroth moments
m0_unreg = [moment_sympy(x, 0) for x in m_unreg]
# Define the first moments
m1_unreg = [moment_sympy(x, 1) for x in m_unreg]
# Define the second moments
m2_unreg = [moment_sympy(x, 2) for x in m_unreg]
# Define the second moments
m3_unreg = [moment_sympy(x, 3) for x in m_unreg]

m0_unreg, m1_unreg, m2_unreg, m3_unreg
```


```python
# define "global" moments
mm0_unreg = moment_sympy('\mathbf{m}_u', 0)
mm1_unreg = moment_sympy('\mathbf{m}_u', None)
mm2_unreg = moment_sympy('\mathbf{m}_u', 2)

mm0_unreg, mm1_unreg, mm2_unreg 
```

Now we will define the state transition matrix $\mathbf{K}^u$, the mRNA production matrix $\mathbf{R}^u_m$ and the mRNA degradation matrix $\mathbf{\Gamma}^u_m$.


```python
# Define the rate constant matrix
Kmat_unreg = sympy.Matrix([[-kp_off, kp_on], 
                           [kp_off, -kp_on]])
# Define the production matrix
Rmat_unreg = sympy.Matrix([[rm, 0], 
                           [0, 0]])
# Define the degradation matrix
Gmat_unreg = sympy.Matrix([[gm, 0],
                           [0, gm]])

Kmat_unreg, Rmat_unreg, Gmat_unreg
```

### Zeroth moment two-state promoter
In order to solve for the two-state promoter we have two equations to use
$$
1 = \mathbf{u} \cdot \left\langle \mathbf{m}^0 \right\rangle,
$$
where $u = (1, 1)$ and $\left\langle \mathbf{m}^0 \right\rangle$ is the vector containing all the states, and
$$
\mathbf{0} = \mathbf{K}^u_m \left\langle \mathbf{m}^0 \right\rangle.
$$

We will use the linear algebra solver where we write an extended matrix with the last column being the right hand side of these equations.

Let's now define this extended matrix.


```python
# Add the extra row for the normalization constraint
K_sys_unreg = Kmat_unreg.row_insert(0, sympy.Matrix([1, 1]).transpose())

# Define right-side array
b = sympy.Matrix([1, 0, 0])

# Define augmented matrix
zero_aug_unreg = K_sys_unreg.col_insert(2, b)
zero_aug_unreg
```

Now we solve the equation for the two state moments.


```python
# Solve linear system for the three zeroth moments
zero_sol_unreg = sympy.solve_linear_system(zero_aug_unreg, *m0_unreg)
zero_sol_unreg
```

### First moment two-state promoter

Now that we have the zeroth moment that indicates the probability of being on either state let's now compute the first moment. For this we have the matrix equation
$$
\mathbf{0} = \left( \mathbf{K}^u_m - \mathbf{\Gamma}^u_m  \right)
             + \mathbf{R}^u_m \left\langle \mathbf{m}^0 \right\rangle.
$$

Let's write this system in the form of the augmented matrix.


```python
# Define augmented matrix for first moment
first_aug_unreg = (Kmat_unreg - Gmat_unreg).col_insert(2, 
                                            -Rmat_unreg * sympy.Matrix(m0_unreg))
first_aug_unreg
```

And again we use the `sympy` solver to solve for the first moments


```python
# Solve linear system for the three moments
first_sol_unreg = sympy.solve_linear_system(first_aug_unreg, *m1_unreg)
first_sol_unreg
```

We now sum both elements and substitute the solution for the zeroth moments


```python
first_unreg = sum(first_sol_unreg.values()).simplify()
first_unreg_s = first_unreg.subs(zero_sol_unreg)
first_unreg, first_unreg_s
```

Let's now define the first moment

### Second moment two-state promoter

For the second moment we have an equation of the form
$$
\mathbf{0} = \left( \mathbf{K}^u_m - 2 \mathbf{\Gamma}^u_n \right) 
             \left\langle \mathbf{m}^2 \right\rangle
             + \left( 2 \mathbf{R}^u_m + \mathbf{\Gamma}^u_m \right)
             \left\langle \mathbf{m}^1 \right\rangle
             + \mathbf{R}^u_m
             \left\langle \mathbf{m}^0 \right\rangle
$$

Let's again define this augmented matrix.


```python
# Define augmented matrix for second moment
second_aug_unreg = (2 * Gmat_unreg - Kmat_unreg).col_insert(2, 
                   (2 * Rmat_unreg + Gmat_unreg) * sympy.Matrix(m1_unreg)\
                   + Rmat_unreg * sympy.Matrix(m0_unreg))
second_aug_unreg
```


```python
# Solve linear system for the three moments
second_sol_unreg = sympy.solve_linear_system(second_aug_unreg, *m2_unreg)
second_sol_unreg
```

Again let's add both of these moments to obtain the global second moment.


```python
second_unreg = sum(second_sol_unreg.values()).simplify()
second_unreg
```

Let's now substitute the solution for the sum of the first moments.


```python
second_unreg_s = second_unreg.subs(gm * m1_unreg[0] + gm * m1_unreg[1], 
                  gm * sympy.simplify(sum(first_sol_unreg.values())))
second_unreg_s.simplify()
```

Let's rewrite this expression in terms of the global moments.


```python
second_unreg_moments = second_unreg.\
                       subs(sum(gm * sympy.Matrix(m1_unreg)), gm * mm1_unreg).\
                       expand().\
                       subs(rm / gm * m0_unreg[0], mm1_unreg)#.\

second_unreg_moments
```

This is the answer shown in PBoC. Now is a matter of substituting the other solutions to obtain the second moment.


```python
second_unreg_ss = sympy.simplify(second_unreg_s).\
                  subs(first_sol_unreg).subs(zero_sol_unreg).\
                  collect(kp_on / (kp_on + kp_off))
second_unreg_ss
```

Let's now compute the variance and the fano-factor to confirm that we obtain the exact same results as derived in PBoC.


```python
# Compute the variance of the unregulated promoter
var_unreg = second_unreg - first_unreg**2
var_unreg_s = var_unreg.subs(first_sol_unreg).subs(zero_sol_unreg)
var_unreg.ratsimp(), var_unreg_s.simplify()
```

Again let's express this as a function of the global moments.


```python
var_unreg_moments = var_unreg.\
                    expand().\
                    subs(sum(m1_unreg) / 2, mm1_unreg / 2).\
                    subs(first_sol_unreg).\
                    subs((rm / gm * m0_unreg[0]), mm1_unreg).\
                    collect(mm1_unreg)

var_unreg_moments
```


```python
# Compute the fano factor of the unregulated promoter
fano_unreg = var_unreg / first_unreg
fano_unreg_s = fano_unreg.subs(first_sol_unreg).subs(zero_sol_unreg)
fano_unreg.ratsimp(), fano_unreg_s.simplify()
```

Again we express this in terms of the global moments


```python
fano_unreg_moments = (var_unreg_moments / mm1_unreg)
fano_unreg_moments
```

Let's manually define the fano-factor as reported in PBoC and compare both expressions.


```python
fano_pboc = 1 + kp_off / (kp_off + kp_on) * rm / (kp_off + kp_on + gm)
fano_pboc
```


```python
fano_pboc.ratsimp() == fano_unreg_s.ratsimp()
```




    True



### Third moment two-state promoter

As an extra derivation we will include the 3rd moment in our calculation. This will be useful to define the equations to compute the parameters. The matrix equation for the 3rd moment is given by
$$
\mathbf{0} = \left( \mathbf{K}^u_m - 3 \mathbf{\Gamma}^u_n \right) 
             \left\langle \mathbf{m}^2 \right\rangle
             + \left( 3 \mathbf{R}^u_m + 3 \mathbf{\Gamma}^u_m \right)
             \left\langle \mathbf{m}^2 \right\rangle
             + \left( 3 \mathbf{R}^u_m - \mathbf{\Gamma}^u_m \right)
             \left\langle \mathbf{m}^1 \right\rangle
             + \mathbf{R}^u_m
             \left\langle \mathbf{m}^0 \right\rangle
$$

Let's again define this augmented matrix.


```python
# Define augmented matrix for third moment
third_aug_unreg = (3 * Gmat_unreg - Kmat_unreg).col_insert(2, 
                  (3 * Rmat_unreg + 3 * Gmat_unreg) * sympy.Matrix(m2_unreg)\
                  + (3 * Rmat_unreg - Gmat_unreg) * sympy.Matrix(m1_unreg)
                  + Rmat_unreg * sympy.Matrix(m0_unreg))
third_aug_unreg
```


```python
# Solve linear system for the three moments
third_sol_unreg = sympy.solve_linear_system(third_aug_unreg, *m3_unreg)
third_sol_unreg
```

As usual let's add these moments


```python
third_unreg = sum(third_sol_unreg.values()).simplify()
third_unreg
```

We again write this in terms of the global momments.


```python
third_unreg_moments = third_unreg.\
                      subs(sum(m1_unreg) / 3, mm1_unreg / 3).\
                      subs(sum(m2_unreg), mm2_unreg).\
                      collect(rm / gm)
#                       subs(first_sol_unreg).\
#                       subs(rm / gm * m0_unreg[0], mm1_unreg).\

third_unreg_moments
```

We also define the third moment with the corresponding substitutions.


```python
third_unreg_ss = third_unreg.subs(second_sol_unreg).\
                 subs(first_sol_unreg).subs(zero_sol_unreg)
third_unreg_ss.simplify()
```

## Deriving expressions for the parameters.

Just as Peccoud & Ycart did we will find equations that solve for the parameters in the equation.

In their paper Peccoud & Ycart note that by computing the exponential moments they can generate a system of equations simple enough to solve for almost all parameters in the equation. Not all four parameters can be uniquely determined since the model is homogeneous in time. What that implies is that multiplying the time scale by a constant is equivalent to dividing all parameters by the same constant and this does not affect the steady-state distribution. But if we define the parameters in units of $\gm$ we can find equations for the other three parameters.

The so-called exponential moments are defined as
$$
e_n \equiv \left\langle m (m - 1) \cdots (m - n + 1) \right\rangle.
$$
Given the linearity property of expected values these moments can be computed from knowing up to the highest regular moment. The first three exponential moments are given by
$$
e_1 = \mm,\\
e_2 = \left\langle m (m - 1) \right\rangle = \left\langle m^2 \right\rangle - \mm,\\
e_3 = \left\langle m (m - 1) (m - 2) \right\rangle =
\left\langle m^3 \right\rangle - 3\left\langle m^2 \right\rangle + 2\mm
$$

Peccoud and Ycart cleverly realize that by defining the successive ratios of these first three exponential moments one can get a system of three equations from which to extract the parameters. The consecutive ratios are simply given by
$$
r_1 \equiv e_1,\\
r_2 \equiv e_2 / e_1,
$$
and
$$
r_3 \equiv e_3 / e_2.
$$

Let's now define these three equations.


```python
r1, r2, r3 = sympy.symbols('r_1 r_2 r_3')
r1_unreg = first_unreg_s.subs(gm, 1)
sympy.Eq(r1, r1_unreg)
```


```python
r2_unreg = ((second_unreg_ss - first_unreg_s) /\
            first_unreg_s).subs(gm, 1).simplify()
sympy.Eq(r2, r2_unreg)
```


```python
r3_unreg = ((third_unreg_ss - 3 * second_unreg_ss + 2 * first_unreg_s) / \
     (second_unreg_ss - first_unreg_s)).subs(gm, 1).simplify()
sympy.Eq(r3, r3_unreg)
```

Now let's use these three equations to solve for the parameters.


```python
infer_param = [kp_on, kp_off, rm]
param_moments = sympy.solve([r1_unreg - r1,
                             r2_unreg - r2,
                             r3_unreg - r3], *infer_param)[0]
param_moments = dict(zip(infer_param, param_moments))
param_moments
```

So by simply computing this experimentally accessible quantities we can obtain estimates for the values of these kinetic parameters.

### Export solutions to `pickle` object

In order to import these results in another notebook we will save the dictionaries with the resutls in a pickle object.


```python
with open(tmpdir + 'two_state_mRNA_solutions.pkl', 'wb') as file:
    pickle.dump(zero_sol_unreg, file)
    pickle.dump(first_sol_unreg, file)
    pickle.dump(second_sol_unreg, file)
    pickle.dump(third_sol_unreg, file)
```


```python
with open(tmpdir + 'two_state_mRNA_moments.pkl', 'wb') as file:
    pickle.dump(first_unreg, file)
    pickle.dump(second_unreg_moments, file)
    pickle.dump(third_unreg_moments, file)
```

### Export `lambdify` functions to `dill` objects

Since we want to keep the notebooks in an organized fashion, the numerical exploration of these functions will be done in a different notebook. Therefore here we will generate `lambdify` functions that will be saved into an object using the `dill` library.

Let' first generate lambdify functions for the moments.


```python
# Define array containing variables
var = [kp_on, kp_off, rm, gm]
# Define lambdify functions for the first three moments
first_unreg_num = sympy.lambdify(var, first_unreg_s)
second_unreg_num = sympy.lambdify(var, second_unreg_ss)
third_unreg_num = sympy.lambdify(var, third_unreg_ss)
```

Now we export the functions into a file using `dill`


```python
# Change dill setting to allow the export of the functions
dill.settings['recurse'] = True

# Open file to save functions
with open(tmpdir + 'two_state_mRNA_lambdify.dill', 'wb') as file:
    dill.dump(first_unreg_num, file)
    dill.dump(second_unreg_num, file)
    dill.dump(third_unreg_num, file)
```

# Three-state promoter (mRNA)

Let's now derive the moments for the three-state promoter. The matrix equations remain the same between the two- and the three-state promoter. The only difference is the definition of the state transition, production and degradation matrix.

Let's first define the moment variables


```python
# Define the string for the states
m_reg = ['m_{P_{reg}}', 'm_{E_{reg}}', 'm_{R_{reg}}']

# Define the zeroth moments
m0_reg = [moment_sympy(x, 0) for x in m_reg]
# Define the first moments
m1_reg = [moment_sympy(x, 1) for x in m_reg]
# Define the second moments
m2_reg = [moment_sympy(x, 2) for x in m_reg]
# Define the second moments
m3_reg = [moment_sympy(x, 3) for x in m_reg]

m0_reg, m1_reg, m2_reg, m3_reg
```


```python
# define "global" moments
mm0_reg = moment_sympy('\mathbf{m}_r', 0)
mm1_reg = moment_sympy('\mathbf{m}_r', None)
mm2_reg = moment_sympy('\mathbf{m}_r', 2)

mm0_reg, mm1_reg, mm2_reg 
```

Now we will define the state transition matrix $\mathbf{K}^r$, the mRNA production matrix $\mathbf{R}^r_m$ and the mRNA degradation matrix $\mathbf{\Gamma}^r_m$.


```python
# Define the rate constant matrix
Kmat_reg = sympy.Matrix([[-kp_off, kp_on, 0], 
                         [kp_off, -(kp_on + kr_on), kr_off],
                         [0, kr_on, -kr_off]])
# Define the production matrix
Rmat_reg = sympy.Matrix([[rm, 0, 0], 
                         [0, 0, 0],
                         [0, 0, 0]])
# Define the degradation matrix
Gmat_reg = sympy.Matrix([[gm, 0, 0],
                         [0, gm, 0],
                         [0, 0, gm]])

Kmat_reg, Rmat_reg, Gmat_reg
```

### Zeroth moment three-state promoter
To solve the zeroth moment for the three-state promoter we use the same set of equations as for the two-state promoter, this is
$$
1 = \mathbf{u} \cdot \left\langle \mathbf{m}^0 \right\rangle,
$$
where this time $u = (1, 1, 1)$ and $\left\langle \mathbf{m}^0 \right\rangle$ is the vector containing all *three* states. We use as well
$$
\mathbf{0} = \mathbf{K}^r_m \left\langle \mathbf{m}^0 \right\rangle.
$$

We will again use the linear algebra solver where we write an extended matrix with the last column being the right hand side of these equations.

Let's now define this extended matrix.


```python
# Add the extra row for the normalization constraint
K_sys_reg = Kmat_reg.row_insert(0, sympy.Matrix([1, 1, 1]).transpose())

# Define right-side array
b = sympy.Matrix([1, 0, 0, 0])

# Define augmented matrix
zero_aug_reg = K_sys_reg.col_insert(3, b)
zero_aug_reg
```

Let's solve the system for our three zeroth moments.


```python
# Solve linear system for the three zeroth moments
zero_sol_reg = sympy.solve_linear_system(zero_aug_reg, *m0_reg)
zero_sol_reg
```

Let's quickly confirm that these moments add up to one since each represents the probability of being on each of the promoter states.


```python
sum(zero_sol_reg.values()).simplify()
```

### First moment three-state promoter

Now that we have the zeroth moment that indicates the probability of being on eeach state let's now compute the first moment. We use the same matrix equation as before, this is
$$
\mathbf{0} = \left( \mathbf{K}^r_m - \mathbf{\Gamma}^r_m  \right)
             + \mathbf{R}^r_m \left\langle \mathbf{m}^0 \right\rangle.
$$

Let's write this system in the form of the augmented matrix.


```python
# Define augmented matrix for first moment
first_aug_reg = (Kmat_reg - Gmat_reg).col_insert(3, 
                                            -Rmat_reg * sympy.Matrix(m0_reg))
first_aug_reg
```

Now let's solve the system.


```python
# Solve linear system for the three moments
first_sol_reg = sympy.solve_linear_system(first_aug_reg, *m1_reg)
first_sol_reg
```

Let's add all these first moments to obtain the complete distribution first moment.


```python
first_reg = sum(first_sol_reg.values()).simplify()
first_reg_s = first_reg.subs(zero_sol_reg)
first_reg, first_reg_s
```

#### fold-change.

A very relevant quantity to compute is the fold-change in gene expression. Since this quantity must be consistent between the thermodynamic framework and the master equation framework this quantity will constraint the value of the parameters that we obtain at the end of the day.


```python
# Compute fold-change
fold_change = first_reg_s.subs(zero_sol_reg) / first_unreg_s
fold_change.factor(kr_off)
```

### Second moment three-state promoter

Finally for the second moment we have again an equation of the form
$$
\mathbf{0} = \left( \mathbf{K}^r_m - 2 \mathbf{\Gamma}^r_n \right) 
             \left\langle \mathbf{m}^2 \right\rangle
             + \left( 2 \mathbf{R}^r_m + \mathbf{\Gamma}^r_m \right)
             \left\langle \mathbf{m}^1 \right\rangle
             + \mathbf{R}^r_m
             \left\langle \mathbf{m}^0 \right\rangle
$$

Let's again define this augmented matrix.


```python
# Define augmented matrix for second moment
second_aug_reg = (2 * Gmat_reg - Kmat_reg).col_insert(3, 
                   (2 * Rmat_reg + Gmat_reg) * sympy.Matrix(m1_reg)\
                   + Rmat_reg * sympy.Matrix(m0_reg))
second_aug_reg
```

And now we solve the second moment.


```python
# Solve linear system for the three moments
second_sol_reg = sympy.solve_linear_system(second_aug_reg, *m2_reg)
second_sol_reg
```

Let's compute the global second moment adding all these second moments.


```python
second_reg = sum(second_sol_reg.values()).simplify().ratsimp()
second_reg_s = second_reg.subs(first_sol_reg).subs(zero_sol_reg).simplify()
second_reg
```

Let's express this in terms of the global moments


```python
second_reg_moments = second_reg.subs(sum(sympy.Matrix(m1_reg) / 2), mm1_reg / 2).\
                                expand().\
                                subs(rm / gm * m0_reg[0], mm1_reg)
second_reg_moments.expand()
```

#### Noise and fano-factor.

Experimentally relevant quantities are the noise (standard deviation / mean) and the fano-factor (variance / mean). Let's compute these quantities for the three-state promoter.

First let's compute the variance.


```python
# Compute the variance for the three-state promoter
var_reg = (second_reg - first_reg**2)
var_reg_s = var_reg.subs(first_sol_reg).subs(zero_sol_reg).simplify()
var_reg.ratsimp()
```

Expressing this in terms of moments gives


```python
var_reg_moments = (second_reg_moments - mm1_reg**2)
var_reg_moments
```

Now the noise squared $\eta_{reg}^2$


```python
# Compute the noise-square for the three-state promoter
noise_square_reg = (var_reg / first_reg**2)
noise_square_reg_s = (var_reg_s / first_reg_s**2).simplify()
noise_square_reg.expand()
```


```python
noise_square_reg_moments = (var_reg_moments / mm1_reg**2).expand()
noise_square_reg_moments
```

And finally the fano-factor.


```python
fano_reg = var_reg / first_reg
fano_reg_s = var_reg_s / first_reg_s
fano_reg.ratsimp()
```


```python
fano_reg_moments = (var_reg_moments / mm1_reg).expand()
fano_reg_moments
```

### Third moment three-state promoter

As an extra derivation we will include the 3rd moment in our calculation. This will be useful to define the equations to compute the parameters. The matrix equation for the 3rd moment is given by
$$
\mathbf{0} = \left( \mathbf{K}^r_m - 3 \mathbf{\Gamma}^r_n \right) 
             \left\langle \mathbf{m}^2 \right\rangle
             + \left( 3 \mathbf{R}^r_m + 3 \mathbf{\Gamma}^r_m \right)
             \left\langle \mathbf{m}^2 \right\rangle
             + \left( 3 \mathbf{R}^r_m - \mathbf{\Gamma}^r_m \right)
             \left\langle \mathbf{m}^1 \right\rangle
             + \mathbf{R}^r_m
             \left\langle \mathbf{m}^0 \right\rangle
$$

Let's again define this augmented matrix.


```python
# Define augmented matrix for third moment
third_aug_reg = (3 * Gmat_reg - Kmat_reg).col_insert(3, 
                  (3 * Rmat_reg + 3 * Gmat_reg) * sympy.Matrix(m2_reg)\
                  + (3 * Rmat_reg - Gmat_reg) * sympy.Matrix(m1_reg)
                  + Rmat_reg * sympy.Matrix(m0_reg))
third_aug_reg
```


```python
# Solve linear system for the three moments
third_sol_reg = sympy.solve_linear_system(third_aug_reg, *m3_reg)
third_sol_reg
```


```python
third_reg = sum(third_sol_reg.values()).simplify()
third_reg_s = third_reg.subs(second_sol_reg).subs(first_sol_reg).\
                        subs(zero_sol_reg).simplify()
third_reg 
```

Once again we express this in terms of the global moments.


```python
third_reg_moments = third_reg.subs(sum(sympy.Matrix(m1_reg) / 3), mm1_reg / 3).\
                              subs(sum(sympy.Matrix(m2_reg)), mm2_reg).\
                              subs(rm / gm * m0_reg[0], mm1_reg).\
                              collect(rm / gm)
third_reg_moments
```

### Export solutions to `pickle` object

In order to import these results in another notebook we will save the dictionaries with the resutls in a pickle object.


```python
with open(tmpdir + 'three_state_mRNA_solutions.pkl', 'wb') as file:
    pickle.dump(first_sol_reg, file)
    pickle.dump(second_sol_reg, file)
    pickle.dump(third_sol_reg, file)
```

### Export `lambdify` functions to `dill` objects

Since we want to keep the notebooks in an organized fashion, the numerical exploration of these functions will be done in a different notebook. Therefore here we will generate `lambdify` functions that will be saved into an object using the `dill` library.

Let' first generate lambdify functions for the moments.


```python
# Define array containing variables
var = [kr_on, kr_off, kp_on, kp_off, rm, gm]
# Define lambdify functions for the first three moments
first_reg_num = sympy.lambdify(var, first_reg_s)
second_reg_num = sympy.lambdify(var, second_reg_s)
third_reg_num = sympy.lambdify(var, third_reg_s)
```

Now let's export them using `dill`.


```python
# Change dill setting to allow the export of the functions
dill.settings['recurse'] = True

# Open file to save functions
with open(tmpdir + 'three_state_mRNA_lambdify.dill', 'wb') as file:
    dill.dump(first_reg_num, file)
    dill.dump(second_reg_num, file)
    dill.dump(third_reg_num, file)
```

### Exploring the use of exponential moments to determine parameters.

Let's explore the possibility of using the exponential moments to determine the missing parameters $\kron$ and $\kroff$. We will use the same strategy as the one Peccoud and Ycart used.

Let's first define the first three exponential moments.


```python
r1 = first_reg.subs(zero_sol_reg)
r1
```


```python
r2 = ((second_reg - first_reg) / first_reg).\
     subs(first_sol_reg).subs(zero_sol_reg).simplify()
r2
```


```python
r3 = ((third_reg - 3 * second_reg + 2 * first_reg) / \
     (second_reg - first_reg)).subs(second_sol_reg).subs(first_sol_reg).\
     subs(zero_sol_reg).simplify()
r3
```


```python
param_reg = sympy.solve([r1 - x1, r2 - x2], kr_off, kr_on)
param_reg
```


```python
test = (param_reg[kr_on] / param_reg[kr_off]).simplify().collect(kp_on / kp_off)
test.subs(x1, r1).simplify()
```


```python
sympy.solve(kr_on / kr_off - test, x1)
```


```python

```

b920f213ebc3b6699cfc2b330e3fc1af3b83b82c
ipynb
Jupyter Notebook
src/theory/sandbox/chemical_master_moments_mRNA.ipynb
RPGroup-PBoC/chann_cap
f2a826166fc2d47c424951c616c46d497ed74b39
MIT
2020-08-21T04:06:12.000Z
2022-02-09T07:36:58.000Z
src/theory/sandbox/chemical_master_moments_mRNA.ipynb
RPGroup-PBoC/chann_cap
f2a826166fc2d47c424951c616c46d497ed74b39
MIT
src/theory/sandbox/chemical_master_moments_mRNA.ipynb
RPGroup-PBoC/chann_cap
f2a826166fc2d47c424951c616c46d497ed74b39
MIT
2020-04-29T17:43:28.000Z
2020-09-09T00:20:16.000Z
Qwen/Qwen-72B
1. YES
2. YES