<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Two-state promoter (protein)1.1&nbsp;&nbsp;Import mRNA moments1.2&nbsp;&nbsp;Zeroth moment two-state promoter1.3&nbsp;&nbsp;First moment two-state promoter (protein)1.4&nbsp;&nbsp;mRNA-protein correlation two-state promoter1.5&nbsp;&nbsp;Second moment two-state promoter1.6&nbsp;&nbsp;Third moment two-state promoter $\ee{\bb{p}^3}$1.6.1&nbsp;&nbsp;$\ee{\bb{m}^2\bb{p}}$1.6.2&nbsp;&nbsp;$\ee{\bb{mp}^2}$1.6.3&nbsp;&nbsp;$\ee{\bb{p}^3}$1.7&nbsp;&nbsp;Export lambdify functions to dill objects2&nbsp;&nbsp;Three-state promoter (protein)2.1&nbsp;&nbsp;Zeroth moment three-state promoter2.2&nbsp;&nbsp;First moment two-state promoter (protein)2.2.1&nbsp;&nbsp;fold-change.2.3&nbsp;&nbsp;mRNA-protein correlation three-state promoter2.4&nbsp;&nbsp;Second moment three-state promoter2.5&nbsp;&nbsp;Third moment three-state promoter $\ee{\bb{p}^3}$2.5.1&nbsp;&nbsp;$\ee{\bb{m}^2\bb{p}}$2.5.2&nbsp;&nbsp;$\ee{\bb{mp}^2}$2.5.3&nbsp;&nbsp;$\ee{\bb{p}^3}$2.6&nbsp;&nbsp;Export lambdify functions to dill objects  The moments of the distribution for protein in a two- and three-state promoter. (c) 2017 <NAME>. 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)<jupyter_code>import pickle
import dill # Library that we will use to export lambdify functions

# 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/'<jupyter_output><empty_output><jupyter_text>$\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}$$\newcommand{ee}[1]{\left\langle 1 \right\rangle}$$\newcommand{bb}[1]{\mathbf{1}}$  Two-state promoter (protein)In this section we will perform the linear algebra to solve for the steady-state moments of the protein 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 protein 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.<jupyter_code># 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')
rp, gp = sympy.symbols('r_p gamma_p')

kp_off, kp_on, kr_off, kr_on, rm, gm, rp, gp<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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}')<jupyter_output><empty_output><jupyter_text>Now let's define the moment variables for the unregulated promoter. We will do it for both the mRNA and the protein since the protein moments depend also on mRNA moments.<jupyter_code># 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
# Define the string for the states
p_unreg = ['p_{P_{unreg}}', 'p_{E_{unreg}}']

# Define the zeroth moments
p0_unreg = [moment_sympy(x, 0) for x in p_unreg]
# Define the first moments
p1_unreg = [moment_sympy(x, 1) for x in p_unreg]
# Define the second moments
p2_unreg = [moment_sympy(x, 2) for x in p_unreg]
# Define the second moments
p3_unreg = [moment_sympy(x, 3) for x in p_unreg]

p0_unreg, p1_unreg, p2_unreg, p3_unreg
# define "global" moments and cross correlations
mm0_unreg = moment_sympy('\mathbf{m}_u', 0)
mm1_unreg = moment_sympy('\mathbf{m}_u', None)
mm2_unreg = moment_sympy('\mathbf{m}_u', 2)
mm3_unreg = moment_sympy('\mathbf{m}_u', 3)
pp0_unreg = moment_sympy('\mathbf{p}_u', 0)
pp1_unreg = moment_sympy('\mathbf{p}_u', None)
pp2_unreg = moment_sympy('\mathbf{p}_u', 2)
pp3_unreg = moment_sympy('\mathbf{p}_u', 3)

# Cross corelations
mmpp_unreg = moment_sympy('\mathbf{mp}_u', None)
mm2pp_unreg = moment_sympy('\mathbf{m}^2\mathbf{p}_u', None)
mmpp2_unreg = moment_sympy('\mathbf{m}\mathbf{p}^2_u', None)

mm0_unreg, mm1_unreg, mm2_unreg, mm3_unreg,\
pp0_unreg, pp1_unreg, pp2_unreg, pp3_unreg,\
mmpp_unreg, mm2pp_unreg, mmpp2_unreg<jupyter_output><empty_output><jupyter_text>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$.<jupyter_code># Define the rate constant matrix
Km_unreg = sympy.Matrix([[-kp_off, kp_on], 
                           [kp_off, -kp_on]])
# Define the mRNA production matrix
Rm_unreg = sympy.Matrix([[rm, 0], 
                         [0, 0]])
# Define the mRNA degradation matrix
Gm_unreg = sympy.Matrix([[gm, 0],
                         [0, gm]])

# Define the protein production matrix
Rp_unreg = sympy.Matrix([[rp, 0], 
                         [0, rp]])
# Define the protein degradation matrix
Gp_unreg = sympy.Matrix([[gp, 0],
                         [0, gp]])


Km_unreg, Rm_unreg, Gm_unreg, Rp_unreg, Gp_unreg<jupyter_output><empty_output><jupyter_text>Import mRNA momentsLet's import the solutions we obtained for the mRNA moments that are saved in a `pickle` object.<jupyter_code>with open(tmpdir + 'two_state_mRNA_solutions.pkl', 'rb') as file:
    zero_sol_unreg_m = pickle.load(file)
    first_sol_unreg_m = pickle.load(file)
    second_sol_unreg_m = pickle.load(file)
    third_sol_unreg_m = pickle.load(file)<jupyter_output><empty_output><jupyter_text>We now import the expressions of the moments in termos of the global moments<jupyter_code>with open(tmpdir + 'two_state_mRNA_moments.pkl', 'rb') as file:
    first_unreg_moments_m = pickle.load(file)
    second_unreg_moments_m = pickle.load(file)
    third_unreg_moments_m = pickle.load(file)<jupyter_output><empty_output><jupyter_text>Zeroth moment two-state promoterIn 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.<jupyter_code># Add the extra row for the normalization constraint
K_sys_unreg = Km_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<jupyter_output><empty_output><jupyter_text>Now we solve the equation for the two state moments.<jupyter_code># Solve linear system for the three zeroth moments
zero_sol_unreg = sympy.solve_linear_system(zero_aug_unreg, *m0_unreg)
zero_sol_unreg<jupyter_output><empty_output><jupyter_text>First moment two-state promoter (protein)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 - \mathbf{\Gamma}^u_p  \right)             + \mathbf{R}^u_p \left\langle \mathbf{m} \right\rangle.$$Let's write this system in the form of the augmented matrix.<jupyter_code># Define augmented matrix for first moment
first_aug_unreg = (Km_unreg - Gp_unreg).col_insert(2, 
                                            -Rp_unreg * sympy.Matrix(m1_unreg))
first_aug_unreg<jupyter_output><empty_output><jupyter_text>And again we use the `sympy` solver to solve for the first moments<jupyter_code># Solve linear system for the first protein moments
first_sol_unreg_p = sympy.solve_linear_system(first_aug_unreg, *p1_unreg)<jupyter_output><empty_output><jupyter_text>We now sum both elements and substitute the solution for the zeroth moments<jupyter_code># Sum the results for the first moments
first_unreg_p = sum(first_sol_unreg_p.values()).simplify()
# Substitute the solutions for the mRNA first and zeroth moment
first_unreg_p_s = first_unreg_p.subs(first_sol_unreg_m).\
                            subs(zero_sol_unreg).simplify()
# Display solution in terms of the mRNA global first moment
first_unreg_p,\
first_unreg_p_s # Disiplay substituted solution.<jupyter_output><empty_output><jupyter_text>Let's display the solution for the first moment in terms of the global moments<jupyter_code>first_unreg_moments_p = first_unreg_p.subs(sum(m1_unreg), mm1_unreg)
sympy.Eq(pp1_unreg, first_unreg_moments_p)<jupyter_output><empty_output><jupyter_text>mRNA-protein correlation two-state promoterThe computation of the second moment $\left\langle \mathbf{p}^2 \right\rangle$ depends on the correlation $\left\langle \mathbf{mp} \right\rangle$. To compute this correlation we have the equation$$\mathbf{0} = \left( \mathbf{K}^u - \mathbf{\Gamma}^u_m -  \mathbf{\Gamma}^u_p  \right)             \left\langle \mathbf{mp} \right\rangle              + \mathbf{R}^u_m \left\langle \mathbf{p} \right\rangle             + \mathbf{R}^u_p \left\langle \mathbf{m}^2 \right\rangle$$Let's define this augmented matrix.<jupyter_code># Define the string for the states
corr_mp_unreg = ['mp_{P_{unreg}}', 'mp_{E_{unreg}}']

# Define the zeroth moments
mp_unreg = [moment_sympy(x, None) for x in corr_mp_unreg]
mp_unreg
# define augmented matrix for correlation
corr_aug_unreg = (Km_unreg - Gm_unreg - Gp_unreg).col_insert(2,
                  - Rm_unreg * sympy.Matrix(p1_unreg)\
                  - Rp_unreg * sympy.Matrix(m2_unreg))
corr_aug_unreg<jupyter_output><empty_output><jupyter_text>Now let's solve the system.<jupyter_code># Solve linear system for the three moments
corr_mp_sol_unreg = sympy.solve_linear_system(corr_aug_unreg, *mp_unreg)<jupyter_output><empty_output><jupyter_text>The sum of these terms gives the correlation we are looking for.<jupyter_code># Sum the terms of the correlation
corr_mp_unreg = sum(corr_mp_sol_unreg.values()).simplify()
# Substitutte relevant solutions
corr_mp_unreg_s = corr_mp_unreg.subs(first_sol_unreg_p).\
                                subs(second_sol_unreg_m).\
                                subs(first_sol_unreg_m).\
                                subs(zero_sol_unreg).simplify()
    
# Display solution without substitutions
corr_mp_unreg<jupyter_output><empty_output><jupyter_text>We now display the solution in terms of the global moments.<jupyter_code># Display correlation as a functiono f the second mRNA moment
corr_moments = corr_mp_unreg.collect(rp).subs(sum(m2_unreg), mm2_unreg)
sympy.Eq(mmpp_unreg, corr_moments)<jupyter_output><empty_output><jupyter_text>LEt's implement the approximation of $\gm \gg \gamma_p$. This would give us<jupyter_code>corr_moments_approx = corr_moments.subs(gm + gp, gm)
corr_moments_approx<jupyter_output><empty_output><jupyter_text>Second moment two-state promoterFor the second moment we have an equation of the form$$\mathbf{0} = \left( \mathbf{K}^u - 2 \mathbf{\Gamma}^u_p \right)              \left\langle \mathbf{p}^2 \right\rangle             + \mathbf{\Gamma}^u_p \left\langle \mathbf{p} \right\rangle             + \mathbf{R}^u_p \left\langle \mathbf{m} \right\rangle             + 2 \mathbf{R}^u_p \left\langle \mathbf{mp} \right\rangle$$Let's again define this augmented matrix.<jupyter_code># Define augmented matrix for second moment
second_aug_unreg = (Km_unreg - 2 * Gp_unreg).col_insert(2, 
                    - Gp_unreg * sympy.Matrix(p1_unreg)\
                    - Rp_unreg * sympy.Matrix(m1_unreg)\
                    - 2 * Rp_unreg * sympy.Matrix(mp_unreg))
second_aug_unreg<jupyter_output><empty_output><jupyter_text>Having defined the matrix let's solve the system.<jupyter_code># Solve linear system for the three moments
second_sol_unreg_p = sympy.solve_linear_system(second_aug_unreg, *p2_unreg)<jupyter_output><empty_output><jupyter_text>Again let's add both of these moments to obtain the global second moment.<jupyter_code># Sum the soltuon of second moments
second_unreg_p = sum(second_sol_unreg_p.values()).simplify()
# Substitute relevant solutions
second_unreg_p_s = second_unreg_p.subs(corr_mp_sol_unreg).\
                                    subs(first_sol_unreg_p).\
                                    subs(second_sol_unreg_m).\
                                    subs(first_sol_unreg_m).\
                                    subs(zero_sol_unreg).simplify()

# Display solution without substitutions
second_unreg_p<jupyter_output><empty_output><jupyter_text>Again let's show the solution using the global moments.The complexity of the substitutions comes from step-by-step manually recognizing patterns.<jupyter_code># Display solution
second_unreg_moments_p = second_unreg_p.\
            collect(gp).collect(rp).\
            subs(sum(p1_unreg), pp1_unreg).\
            subs(sum(m1_unreg), mm1_unreg).\
            subs(sum(2 * mp_unreg), 2 * mmpp_unreg).\
            expand().\
            subs(rp / (2 * gp) * mm1_unreg, pp1_unreg / 2)
sympy.Eq(pp2_unreg,  second_unreg_moments_p)<jupyter_output><empty_output><jupyter_text>Using the solution for $\left\langle \mathbf{mp} \right\rangle$ with the $\gm \gg \gamma_p$ approximation gives<jupyter_code>second_unreg_moments_approx = second_unreg_moments_p.\
                              subs(mmpp_unreg, corr_moments_approx).\
                              collect(rp / (gm * gp))

sympy.Eq(pp2_unreg, second_unreg_moments_approx)<jupyter_output><empty_output><jupyter_text>Third moment two-state promoter $\ee{\bb{p}^3}$  $\ee{\bb{m}^2\bb{p}}$To compute the third moment $\ee{\bb{p}^3}$ we need to compute the correlation term $\ee{\bb{mp}^2}$. This correlation needs at the same time the correlation $\ee{\bb{m}^2\bb{p}}$. So we will solve for these two terms.For the correlation $\ee{\bb{m}^2\bb{p}}$ we have the equation$$\bb{0} = \left( \bb{K}^u - 2 \bb{\Gamma}_p^u - \bb{\Gamma}_p^u \right) \ee{\bb{m}^2\bb{p}} +\bb{R}_m^u \ee{\bb{p}} +\left( 2 \bb{R}_m^u + \bb{\Gamma}_m^u \right) \ee{\bb{mp}} +\bb{R}_p^u \ee{\bb{m}^3}.$$Let's define this augmented matrix.<jupyter_code># Define the string for the states
corr_mp_unreg = ['m^2p_{P_{unreg}}', 'm^2p_{E_{unreg}}']

# Define the zeroth moments
m2p_unreg = [moment_sympy(x, None) for x in corr_mp_unreg]
m2p_unreg
# define augmented matrix for correlation
m2p_aug_unreg = (Km_unreg - 2 * Gm_unreg - Gp_unreg).col_insert(2,
                  - Rm_unreg * sympy.Matrix(p1_unreg)\
                  - (2 * Rm_unreg + Gm_unreg) * sympy.Matrix(mp_unreg)\
                  - Rp_unreg * sympy.Matrix(m3_unreg))
m2p_aug_unreg<jupyter_output><empty_output><jupyter_text>Let's now solve the system.<jupyter_code># Solve linear system for the three moments
m2p_sol_unreg = sympy.solve_linear_system(m2p_aug_unreg, *m2p_unreg)<jupyter_output><empty_output><jupyter_text>By adding the terms we can find the correlation we are looking for<jupyter_code># Sum the terms of the m2pelation
corr_m2p_unreg = sum(m2p_sol_unreg.values()).simplify()
# Substitutte relevant solutions
corr_m2p_unreg_s = corr_m2p_unreg.subs(corr_mp_sol_unreg).\
                                    subs(first_sol_unreg_p).\
                                    subs(third_sol_unreg_m).\
                                    subs(second_sol_unreg_m).\
                                    subs(first_sol_unreg_m).\
                                    subs(zero_sol_unreg)
    
# Display solution without substitutions
corr_m2p_unreg<jupyter_output><empty_output><jupyter_text>Let's now display this solution as a function of the global moments.<jupyter_code># Display m2pelation as a functiono f the second mRNA moment
m2p_moments = corr_m2p_unreg.collect(gm).subs(sum(mp_unreg), mmpp_unreg)\
                       .collect(rp).subs(sum(m3_unreg), mm3_unreg)\
                       .collect(rm)
sympy.Eq(mm2pp_unreg, m2p_moments)<jupyter_output><empty_output><jupyter_text>$\ee{\bb{mp}^2}$We have a solution for all these moments, so we can procced to compute the correlation $\ee{\bb{mp}^2}$. For that we have the equation$$\bb{0} = \left( \bb{K}^u - \bb{\Gamma}_m^u - 2 \bb{\Gamma}_p^u \right)\ee{\bb{mp}^2} +\bb{R}_m^u \ee{\bb{p}^2} +2 \bb{R}_p^u \ee{\bb{m}^2\bb{p}} +\bb{R}_p^u \ee{\bb{m}^2} +\bb{\Gamma}_p^u \ee{\bb{mp}}.$$Let's define the agumented matrix for this correlation term.<jupyter_code># Define the string for the states
corr_mp_unreg = ['mp^2_{P_{unreg}}', 'mp^2_{E_{unreg}}']

# Define the zeroth moments
mp2_unreg = [moment_sympy(x, None) for x in corr_mp_unreg]
mp2_unreg
# define augmented matrix for correlation
mp2_aug_unreg = (Km_unreg - Gm_unreg - 2 * Gp_unreg).col_insert(2,
                  - Rm_unreg * sympy.Matrix(p2_unreg)\
                  - 2 * Rp_unreg * sympy.Matrix(m2p_unreg)\
                  - Rp_unreg * sympy.Matrix(m2_unreg)\
                  - Gp_unreg * sympy.Matrix(mp_unreg))
mp2_aug_unreg<jupyter_output><empty_output><jupyter_text>Let's solve the system.<jupyter_code># Solve linear system for the three moments
mp2_sol_unreg = sympy.solve_linear_system(mp2_aug_unreg, *mp2_unreg)<jupyter_output><empty_output><jupyter_text>By adding the terms we can find the correlation we are looking for<jupyter_code># Sum the terms of the mp2elation
corr_mp2_unreg = sum(mp2_sol_unreg.values()).simplify()
# Substitutte relevant solutions
corr_mp2_unreg_s = corr_mp2_unreg.subs(m2p_sol_unreg).\
                                    subs(second_sol_unreg_p).\
                                    subs(corr_mp_sol_unreg).\
                                    subs(first_sol_unreg_p).\
                                    subs(third_sol_unreg_m).\
                                    subs(second_sol_unreg_m).\
                                    subs(first_sol_unreg_m).\
                                    subs(zero_sol_unreg).simplify()
    
# Display solution without substitutions
corr_mp2_unreg<jupyter_output><empty_output><jupyter_text>Let's now display this solution as a function of the global moments.<jupyter_code># Display correlation as a function of the second mRNA moment
mp2_moments = corr_mp2_unreg.collect(gp).subs(sum(mp_unreg), mmpp_unreg)\
                       .collect(rp).subs(sum(m2_unreg), mm2_unreg)\
                       .subs(sum(2 * m2p_unreg), 2 * mm2pp_unreg)
sympy.Eq(mmpp2_unreg, mp2_moments)<jupyter_output><empty_output><jupyter_text>$\ee{\bb{p}^3}$Finally having these two correlations in hand we can compute the third moment of the protein distribution $\ee{\mm{p}^3}$. For that we have the equation$$\bb{0} = \left( \bb{K}^u - 3 \bb{\Gamma}_p^u \right) \ee{\bb{p}^3} +3 \bb{\Gamma}_p^u \ee{\bb{p}^2} -\bb{\Gamma}_p^u \ee{\bb{p}} +3 \bb{R}_p^u \ee{\bb{mp}^2} +3 \bb{R}_p^u \ee{\bb{mp}} +\bb{R}_p^u \ee{\bb{m}}.$$Let's define this agumented matrix.<jupyter_code># define augmented matrix for correlation
p3_aug_unreg = (Km_unreg - 3 * Gp_unreg).col_insert(2,
                  - 3 * Gp_unreg * sympy.Matrix(p2_unreg)\
                  - Gp_unreg * sympy.Matrix(p1_unreg)\
                  - 3 * Rp_unreg * sympy.Matrix(mp2_unreg)\
                  - 3 * Rp_unreg * sympy.Matrix(mp_unreg)\
                  - Rp_unreg * sympy.Matrix(m1_unreg))
p3_aug_unreg<jupyter_output><empty_output><jupyter_text>Let's solve the linear system<jupyter_code>third_sol_unreg_p = sympy.solve_linear_system(p3_aug_unreg, *p3_unreg)<jupyter_output><empty_output><jupyter_text>Now let's add the solutions to obtain the desired moment.<jupyter_code># Sum the soltuon of third moments
third_unreg_p = sum(third_sol_unreg_p.values()).simplify()
# Substitute relevant solutions
third_unreg_p_s = third_unreg_p.subs(mp2_sol_unreg).\
                    subs(m2p_sol_unreg).\
                    subs(second_sol_unreg_p).\
                    subs(corr_mp_sol_unreg).\
                    subs(first_sol_unreg_p).\
                    subs(third_sol_unreg_m).\
                    subs(second_sol_unreg_m).\
                    subs(first_sol_unreg_m).\
                    subs(zero_sol_unreg)

# Display solution without substitutions
third_unreg_p<jupyter_output><empty_output><jupyter_text>Let's write down this in terms of the global moments.<jupyter_code># Display solution
third_unreg_moments_p = third_unreg_p.\
            collect(gp).\
            subs(sum(p1_unreg) / 3, pp1_unreg / 3).\
            subs(sum(p2_unreg), pp2_unreg).\
            collect(rp).\
            subs(sum(m1_unreg) / 3, mm1_unreg / 3).\
            subs(sum(mp2_unreg), mmpp2_unreg).\
            subs(sum(mp_unreg), mmpp_unreg)
sympy.Eq(pp2_unreg,  third_unreg_moments_p)<jupyter_output><empty_output><jupyter_text>--- Export `lambdify` functions to `dill` objectsSince 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.<jupyter_code># Define array containing variables
var = [kp_on, kp_off, rm, gm, rp, gp]
# Define lambdify functions for the first three moments
first_unreg_num = sympy.lambdify(var, first_unreg_p_s)
second_unreg_num = sympy.lambdify(var, second_unreg_p_s)
third_unreg_num = sympy.lambdify(var, third_unreg_p_s)

# Define lambdify functions for the cross correlation terms
corr_mp_unreg_num = sympy.lambdify(var, corr_mp_unreg_s)
corr_m2p_unreg_num = sympy.lambdify(var, corr_m2p_unreg_s)
corr_mp2_unreg_num = sympy.lambdify(var, corr_mp2_unreg_s)
# Change dill setting to allow the export of the functions
dill.settings['recurse'] = True

# Open file to save functions
with open(tmpdir + 'two_state_protein_lambdify.dill', 'wb') as file:
    dill.dump(first_unreg_num, file)
    dill.dump(second_unreg_num, file)
    dill.dump(third_unreg_num, file)
    dill.dump(corr_mp_unreg_num, file)
    dill.dump(corr_m2p_unreg_num, file)
    dill.dump(corr_mp2_unreg_num, file)<jupyter_output><empty_output><jupyter_text>--- Three-state promoter (protein)--- 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<jupyter_code># 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
# Define the string for the states
p_reg = ['p_{P_{reg}}', 'p_{E_{reg}}', 'p_{R_{reg}}']

# Define the zeroth moments
p0_reg = [moment_sympy(x, 0) for x in p_reg]
# Define the first moments
p1_reg = [moment_sympy(x, 1) for x in p_reg]
# Define the second moments
p2_reg = [moment_sympy(x, 2) for x in p_reg]
# Define the third moments
p3_reg = [moment_sympy(x, 3) for x in p_reg]

p0_reg, p1_reg, p2_reg, p3_reg
# 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)
mm3_reg = moment_sympy('\mathbf{m}_r', 3)
pp0_reg = moment_sympy('\mathbf{p}_r', 0)
pp1_reg = moment_sympy('\mathbf{p}_r', None)
pp2_reg = moment_sympy('\mathbf{p}_r', 2)
pp3_reg = moment_sympy('\mathbf{p}_r', 3)

# Cross corelations
mmpp_reg = moment_sympy('\mathbf{mp}_r', None)
mm2pp_reg = moment_sympy('\mathbf{m}^2\mathbf{p}_r', None)
mmpp2_reg = moment_sympy('\mathbf{m}\mathbf{p}^2_r', None)

mm0_reg, mm1_reg, mm2_reg, mm3_reg,\
pp0_reg, pp1_reg, pp2_reg, pp3_reg,\
mmpp_reg, mm2pp_reg, mmpp2_reg
with open(tmpdir + 'three_state_mRNA_solutions.pkl', 'rb') as file:
    first_sol_reg_m = pickle.load(file)
    second_sol_reg_m = pickle.load(file)
    third_sol_reg_m = pickle.load(file)
# Define the rate constant matrix
Km_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
Rm_reg = sympy.Matrix([[rm, 0, 0],
                       [0, 0, 0],
                       [0, 0, 0]])
# Define the degradation matrix
Gm_reg = sympy.Matrix([[gm, 0, 0],
                       [0, gm, 0],
                       [0, 0, gm]])

# Define the production matrix
Rp_reg = sympy.Matrix([[rp, 0, 0],
                       [0, rp, 0],
                       [0, 0, rp]])

# Define the production matrix
Gp_reg = sympy.Matrix([[gp, 0, 0],
                       [0, gp, 0],
                       [0, 0, gp]])


Km_reg, Rm_reg, Gm_reg, Rp_reg, Gp_reg<jupyter_output><empty_output><jupyter_text>Zeroth moment three-state promoterTo 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 \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.<jupyter_code># Add the extra row for the normalization constraint
K_sys_reg = Km_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<jupyter_output><empty_output><jupyter_text>Let's solve the system for our three zeroth moments.<jupyter_code># Solve linear system for the three zeroth moments
zero_sol_reg = sympy.solve_linear_system(zero_aug_reg, *m0_reg)
zero_sol_reg<jupyter_output><empty_output><jupyter_text>Let's quickly confirm that these moments add up to one since each represents the probability of being on each of the promoter states.<jupyter_code>sum(zero_sol_reg.values()).simplify()<jupyter_output><empty_output><jupyter_text>First moment two-state promoter (protein)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}^r - \mathbf{\Gamma}^r_p  \right)             + \mathbf{R}^r_p \left\langle \mathbf{m} \right\rangle.$$Let's write this system in the form of the augmented matrix.<jupyter_code># Define augmented matrix for first moment
first_aug_reg = (Km_reg - Gp_reg).col_insert(3, 
                -Rp_reg * sympy.Matrix(m1_reg))
first_aug_reg<jupyter_output><empty_output><jupyter_text>Now let's solve the system.<jupyter_code># Solve linear system for the first protein moments
first_sol_reg_p = sympy.solve_linear_system(first_aug_reg, *p1_reg)<jupyter_output><empty_output><jupyter_text>Let's sum the elements of this solution to compute the global moment.<jupyter_code># Sum the results for the first moments
first_reg_p = sum(first_sol_reg_p.values()).simplify()
# Substitute the solutions for the mRNA first and zeroth moment
first_reg_p_s = first_reg_p.subs(first_sol_reg_m).\
                            subs(zero_sol_reg).simplify()
# Display solution in terms of the mRNA global first moment
first_reg_p,\
first_reg_p_s # Disiplay substituted solution.<jupyter_output><empty_output><jupyter_text>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.<jupyter_code># Compute fold-change
fold_change = first_reg_p_s / first_unreg_p_s
fold_change.factor(kr_off)<jupyter_output><empty_output><jupyter_text>This is as expected the same expression that we obtained at the mRNA level.  mRNA-protein correlation three-state promoterThe computation of the second moment $\left\langle \mathbf{p}^2 \right\rangle$ depends on the correlation $\left\langle \mathbf{mp} \right\rangle$. To compute this correlation we have the equation$$\mathbf{0} = \left( \mathbf{K}^r - \mathbf{\Gamma}^r_m -  \mathbf{\Gamma}^r_p  \right)             \left\langle \mathbf{mp} \right\rangle              + \mathbf{R}^r_m \left\langle \mathbf{p} \right\rangle             + \mathbf{R}^r_p \left\langle \mathbf{m}^2 \right\rangle$$Let's define this augmented matrix.<jupyter_code># Define the string for the states
corr_mp_reg = ['mp_{P_{reg}}', 'mp_{E_{reg}}', 'mp_{R_{reg}}']

# Define the zeroth moments
mp_reg = [moment_sympy(x, None) for x in corr_mp_reg]
mp_reg
# define augmented matrix for correlation
corr_aug_reg = (Km_reg - Gm_reg - Gp_reg).col_insert(3,
                - Rm_reg * sympy.Matrix(p1_reg)\
                - Rp_reg * sympy.Matrix(m2_reg))
corr_aug_reg
# Solve linear system for the three moments
corr_mp_sol_reg = sympy.solve_linear_system(corr_aug_reg, *mp_reg)<jupyter_output><empty_output><jupyter_text>Now let's add each of these solutions to obtain the total correlation.<jupyter_code># Sum the terms of the correlation
corr_mp_reg = sum(corr_mp_sol_reg.values()).simplify()
# Substitutte relevant solutions
corr_mp_reg_s = corr_mp_reg.subs(first_sol_reg_p).subs(second_sol_reg_m).\
               subs(first_sol_reg_m).subs(zero_sol_reg)
    
# Display solution without substitutions
corr_mp_reg<jupyter_output><empty_output><jupyter_text>Let's now write this in terms of the global moments<jupyter_code># Display correlation as a functiono f the second mRNA moment
corr_moments = corr_mp_reg.collect(rp).subs(sum(m2_reg), mm2_reg)
sympy.Eq(mmpp_reg, corr_moments)<jupyter_output><empty_output><jupyter_text>Second moment three-state promoterFor the second moment we have an equation of the form$$\mathbf{0} = \left( \mathbf{K}^r - 2 \mathbf{\Gamma}^r_p \right)              \left\langle \mathbf{p}^2 \right\rangle             + \mathbf{\Gamma}^r_p \left\langle \mathbf{p} \right\rangle             + \mathbf{R}^r_p \mm             + 2 \mathbf{R}^r_p \left\langle \mathbf{mp} \right\rangle$$Let's again define this augmented matrix.<jupyter_code># Define augmented matrix for second moment
second_aug_reg = (Km_reg - 2 * Gp_reg).col_insert(3, 
                    - Gp_reg * sympy.Matrix(p1_reg)\
                    - Rp_reg * sympy.Matrix(m1_reg)\
                    - 2 * Rp_reg * sympy.Matrix(mp_reg))
second_aug_reg<jupyter_output><empty_output><jupyter_text>Having defined the augmented matrix let's solve the linear system.<jupyter_code># Solve linear system for the three moments
second_sol_reg_p = sympy.solve_linear_system(second_aug_reg, *p2_reg)<jupyter_output><empty_output><jupyter_text>Now we add all these elements to obtain the global second moment<jupyter_code># Sum the soltuon of second moments
second_reg_p = sum(second_sol_reg_p.values()).simplify()
# Substitute relevant solutions
second_reg_p_s = second_reg_p.subs(corr_mp_sol_reg).subs(first_sol_reg_p).\
                 subs(second_sol_reg_m).subs(first_sol_reg_m).\
                 subs(zero_sol_reg)

# Display solution without substitutions
second_reg_p<jupyter_output><empty_output><jupyter_text>Let's now express this moment as a funciton of the global moments.<jupyter_code># Display solution
second_reg_moments_p = second_reg_p.\
            collect(gp).collect(rp).\
            subs(sum(p1_reg), pp1_reg).\
            subs(sum(m1_reg), mm1_reg).\
            subs(sum(2 * mp_reg), 2 * mmpp_reg).\
            expand().\
            subs(rp / (2 * gp) * mm1_reg, pp1_reg / 2)
sympy.Eq(pp2_reg,  second_reg_moments_p)<jupyter_output><empty_output><jupyter_text>It has the exact same funcitonal form as the two-state promoter. The only difference is obviously the definition of each of the moments.  Third moment three-state promoter $\ee{\bb{p}^3}$ $\ee{\bb{m}^2\bb{p}}$To compute the third moment $\ee{\bb{p}^3}$ we need to compute the correlation term $\ee{\bb{mp}^2}$. This correlation needs at the same time the correlation $\ee{\bb{m}^2\bb{p}}$. So we will solve for these two terms.For the correlation $\ee{\bb{m}^2\bb{p}}$ we have the equation$$\bb{0} = \left( \bb{K}^r - 2 \bb{\Gamma}_p^r - \bb{\Gamma}_p^r \right) \ee{\bb{m}^2\bb{p}} +\bb{R}_m^r \ee{\bb{p}} +\left( 2 \bb{R}_m^r + \bb{\Gamma}_m^r \right) \ee{\bb{mp}} +\bb{R}_p^r \ee{\bb{m}^3}.$$Let's define this augmented matrix.<jupyter_code># Define the string for the states
corr_mp_reg = ['m^2p_{P_{reg}}', 'm^2p_{E_{reg}}', 'm^2p_{R_{reg}}']

# Define the zeroth moments
m2p_reg = [moment_sympy(x, None) for x in corr_mp_reg]
m2p_reg
# define augmented matrix for correlation
m2p_aug_reg = (Km_reg - 2 * Gm_reg - Gp_reg).col_insert(3,
                  - Rm_reg * sympy.Matrix(p1_reg)\
                  - (2 * Rm_reg + Gm_reg) * sympy.Matrix(mp_reg)\
                  - Rp_reg * sympy.Matrix(m3_reg))
m2p_aug_reg<jupyter_output><empty_output><jupyter_text>Let's now solve the system.<jupyter_code># Solve linear system for the three moments
m2p_sol_reg = sympy.solve_linear_system(m2p_aug_reg, *m2p_reg)<jupyter_output><empty_output><jupyter_text>Let's add these terms to see the total correlation.<jupyter_code># Sum the terms of the m2pelation
corr_m2p_reg = sum(m2p_sol_reg.values()).simplify()
# Substitutte relevant solutions
corr_m2p_reg_s = corr_m2p_reg.subs(second_sol_reg_p).\
                    subs(corr_mp_sol_reg).\
                    subs(first_sol_reg_p).\
                    subs(third_sol_reg_m).\
                    subs(second_sol_reg_m).\
                    subs(first_sol_reg_m).\
                    subs(zero_sol_reg)
    
# Display solution without substitutions
corr_m2p_reg<jupyter_output><empty_output><jupyter_text>Now let's express this in terms of the global moments<jupyter_code># Display m2pelation as a functiono f the second mRNA moment
m2p_moments_reg = corr_m2p_reg.collect(gm).subs(sum(mp_reg), mmpp_reg)\
                       .collect(rp).subs(sum(m3_reg), mm3_reg)\
                       .collect(rm)
sympy.Eq(mm2pp_reg, m2p_moments_reg)<jupyter_output><empty_output><jupyter_text>$\ee{\bb{mp}^2}$We have a solution for all these moments, so we can procced to compute the correlation $\ee{\bb{mp}^2}$. For that we have the equation$$\bb{0} = \left( \bb{K}^r - \bb{\Gamma}_m^r - 2 \bb{\Gamma}_p^r \right)\ee{\bb{mp}^2} +\bb{R}_m^r \ee{\bb{p}^2} +2 \bb{R}_p^r \ee{\bb{m}^2\bb{p}} +\bb{R}_p^r \ee{\bb{m}^2} +\bb{\Gamma}_p^r \ee{\bb{mp}}.$$Let's define the agumented matrix for this correlation term.<jupyter_code># Define the string for the states
corr_mp_reg = ['mp^2_{P_{reg}}', 'mp^2_{E_{reg}}', 'mp^2_{R_{reg}}']

# Define the zeroth moments
mp2_reg = [moment_sympy(x, None) for x in corr_mp_reg]
mp2_reg
# define augmented matrix for correlation
mp2_aug_reg = (Km_reg - Gm_reg - 2 * Gp_reg).col_insert(3,
                  - Rm_reg * sympy.Matrix(p2_reg)\
                  - 2 * Rp_reg * sympy.Matrix(m2p_reg)\
                  - Rp_reg * sympy.Matrix(m2_reg)\
                  - Gp_reg * sympy.Matrix(mp_reg))
mp2_aug_reg<jupyter_output><empty_output><jupyter_text>Now we can solve the system<jupyter_code># Solve linear system for the three moments
mp2_sol_reg = sympy.solve_linear_system(mp2_aug_reg, *mp2_reg)<jupyter_output><empty_output><jupyter_text>Again we add the individual terms to compute the total correlation term.<jupyter_code># Sum the terms of the mp2elation
corr_mp2_reg = sum(mp2_sol_reg.values()).simplify()
# Substitutte relevant solutions
corr_mp2_reg_s = corr_mp2_reg.subs(m2p_sol_reg).\
                    subs(second_sol_reg_p).\
                    subs(corr_mp_sol_reg).\
                    subs(first_sol_reg_p).\
                    subs(third_sol_reg_m).\
                    subs(second_sol_reg_m).\
                    subs(first_sol_reg_m).\
                    subs(zero_sol_reg)
    
# Display solution without substitutions
corr_mp2_reg<jupyter_output><empty_output><jupyter_text>In terms of the total correlation terms this gives<jupyter_code># Display correlation as a function of the second mRNA moment
mp2_moments_reg = corr_mp2_reg.collect(gp).subs(sum(mp_reg), mmpp_reg)\
                       .collect(rp).subs(sum(m2_reg), mm2_reg)\
                       .subs(sum(2 * m2p_reg), 2 * mm2pp_reg)
sympy.Eq(mmpp2_reg, mp2_moments_reg)<jupyter_output><empty_output><jupyter_text>$\ee{\bb{p}^3}$Finally having these two correlations in hand we can compute the third moment of the protein distribution $\ee{\mm{p}^3}$. For that we have the equation$$\bb{0} = \left( \bb{K}^r - 3 \bb{\Gamma}_p^r \right) \ee{\bb{p}^3} +3 \bb{\Gamma}_p^r \ee{\bb{p}^2} -\bb{\Gamma}_p^r \ee{\bb{p}} +3 \bb{R}_p^r \ee{\bb{mp}^2} +3 \bb{R}_p^r \ee{\bb{mp}} +\bb{R}_p^r \ee{\bb{m}}.$$Let's define this agumented matrix.<jupyter_code># define augmented matrix for correlation
p3_aug_reg = (Km_reg - 3 * Gp_reg).col_insert(3,
                  - 3 * Gp_reg * sympy.Matrix(p2_reg)\
                  - Gp_reg * sympy.Matrix(p1_reg)\
                  - 3 * Rp_reg * sympy.Matrix(mp2_reg)\
                  - 3 * Rp_reg * sympy.Matrix(mp_reg)\
                  - Rp_reg * sympy.Matrix(m1_reg))
p3_aug_reg<jupyter_output><empty_output><jupyter_text>We now solve the system.<jupyter_code>third_sol_reg_p = sympy.solve_linear_system(p3_aug_reg, *p3_reg)<jupyter_output><empty_output><jupyter_text>And compute the total third moment again by adding the individual terms.<jupyter_code># Sum the soltuon of third moments
third_reg_p = sum(third_sol_reg_p.values()).simplify()
# Substitute relevant solutions
third_reg_p_s = third_reg_p.subs(mp2_sol_reg).\
                    subs(m2p_sol_reg).\
                    subs(second_sol_reg_p).\
                    subs(corr_mp_sol_reg).\
                    subs(first_sol_reg_p).\
                    subs(third_sol_reg_m).\
                    subs(second_sol_reg_m).\
                    subs(first_sol_reg_m).\
                    subs(zero_sol_reg)

# Display solution without substitutions
third_reg_p<jupyter_output><empty_output><jupyter_text>In term of the total moments this is given by<jupyter_code># Display solution
third_reg_moments_p = third_reg_p.\
            collect(gp).\
            subs(sum(p1_reg) / 3, pp1_reg / 3).\
            subs(sum(p2_reg), pp2_reg).\
            collect(rp).\
            subs(sum(m1_reg) / 3, mm1_reg / 3).\
            subs(sum(mp2_reg), mmpp2_reg).\
            subs(sum(mp_reg), mmpp_reg)
sympy.Eq(pp2_reg,  third_reg_moments_p)<jupyter_output><empty_output><jupyter_text>Export `lambdify` functions to `dill` objectsSince 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.<jupyter_code># Define array containing variables
var = [kr_on, kr_off, kp_on, kp_off, rm, gm, rp, gp]
# Define lambdify functions for the first three moments
first_reg_num = sympy.lambdify(var, first_reg_p_s)
second_reg_num = sympy.lambdify(var, second_reg_p_s)
third_reg_num = sympy.lambdify(var, third_reg_p_s)

# Define lambdify functions for the cross correlation terms
corr_mp_reg_num = sympy.lambdify(var, corr_mp_reg_s)
corr_m2p_reg_num = sympy.lambdify(var, corr_m2p_reg_s)
corr_mp2_reg_num = sympy.lambdify(var, corr_mp2_reg_s)
# Change dill setting to allow the export of the functions
dill.settings['recurse'] = True

# Open file to save functions
with open(tmpdir + 'three_state_protein_lambdify.dill', 'wb') as file:
    dill.dump(first_reg_num, file)
    dill.dump(second_reg_num, file)
    dill.dump(third_reg_num, file)
    dill.dump(corr_mp_reg_num, file)
    dill.dump(corr_m2p_reg_num, file)
    dill.dump(corr_mp2_reg_num, file)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11791968274445161, 'fraction_numerical': 0.013402135741417616, 'mean_word_length': 3.8404255319148937, 'pattern_counts': {'":': 0, '<': 260, '<?xml version=': 0, '>': 260, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '64', 'dir': 'jupyter-structured-clean-dedup', 'id': '71081', 'license': 'MIT', 'n_tokens_mistral': 15383, 'n_tokens_neox': 14711, 'n_words': 3519, 'path': 'src/theory/sandbox/chemical_master_moments_protein.ipynb', 'repo_name': 'RPGroup-PBoC/chann_cap', 'size': '41864'}