```python
from sympy import *
import numpy as np
from numpy import linalg as LA

import pandas as pd
from pandas import DataFrame, Series
import seaborn as sns
import matplotlib.pyplot as plt
import scipy.optimize as opt
```


```python
# This function calculates the Maxium Lyapunov Exponent for a given steady state

def calculate_MLE(steady_state, n_val, m_val, kc_val, kt_val, gamma_val, gamma_R_val, k_RS_val, beta_N_val,  beta_L_val, beta_R_val):
    
    # First we initialize the parameters and variables in the equations
    # N: Notch receptor
    # L: Ligand (jag1a)
    # R: Repressor (her6/her9)
    # n after N, L or R refers to 'neighbors'
    # gamma are degradation rates
    # beta are production rates
    # kc and ktrans correspond to interaction strength in cis and trans
    # n and m are the hill coefficients (n = m = 1 means no cooperativity)
    N, L, R, Nn, Ln, Rn, n, m, kc, kt, gamma, gamma_R, k_RS, beta_N,  beta_L, beta_R = symbols('N L R Nn Ln Rn n m kc kt gamma gamma_R k_RS beta_N beta_L beta_R')#, real = True)

    N_homo = steady_state[0] # Notch amount at the homogeneous steady state
    L_homo = steady_state[1] # Ligand amount at the homogeneous steady state
    R_homo = steady_state[2] # Repressor amount at the homogeneous steady state
    
    # We define the three equations for dN/dt, dL/dt and dR/dt
    eq1 = beta_N - gamma*N - kt*N*Ln - kc*N*L
    eq2 = beta_L*(1/(1+R**m)) - gamma*L - kt*Nn*L - kc*N*L
    eq3 = beta_R*((N*Ln)**n)/(k_RS + (N*Ln)**n) - gamma_R*R

    # We substitute the parameter values
    eq1_subs = eq1.subs([(beta_N,beta_N_val),(gamma,gamma_val),(kt,kt_val),(kc,kc_val)])
    eq2_subs = eq2.subs([(beta_L,beta_L_val),(gamma,gamma_val),(kt,kt_val),(kc,kc_val), (m, m_val)])
    eq3_subs = eq3.subs([(beta_R,beta_R_val),(gamma_R,gamma_R_val),(k_RS, k_RS_val), (n, n_val)])
    
    #We differentiate the three equations with respect to N, L and R    
    a1 = diff(eq1_subs, N).subs([(N,N_homo), (L,L_homo),             (Ln,L_homo)])
    b1 = diff(eq1_subs, L).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Ln,L_homo)])
    c1 = diff(eq1_subs, R).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Ln,L_homo)])

    a2 = diff(eq2_subs, N).subs([(N,N_homo), (L,L_homo),             (Nn,N_homo)])
    b2 = diff(eq2_subs, L).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Nn,N_homo)])
    c2 = diff(eq2_subs, R).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Nn,N_homo)])

    a3 = diff(eq3_subs, N).subs([(N,N_homo), (L,L_homo),             (Ln,L_homo)])
    b3 = diff(eq3_subs, L).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Ln,L_homo)])
    c3 = diff(eq3_subs, R).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Ln,L_homo)])
    
    # And store these values into the 'H' array
    H = np.array([[a1, b1, c1], [a2, b2, c2], [a3, b3, c3]], dtype = 'float64')

    # Then we differentiate with respect to Nn, Ln and Rn
    a1 = diff(eq1_subs, Nn).subs([(N,N_homo), (L,L_homo),             (Ln,L_homo)])
    b1 = diff(eq1_subs, Ln).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Ln,L_homo)])
    c1 = diff(eq1_subs, Rn).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Ln,L_homo)])

    a2 = diff(eq2_subs, Nn).subs([(N,N_homo), (L,L_homo),             (Nn,N_homo)])
    b2 = diff(eq2_subs, Ln).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Nn,N_homo)])
    c2 = diff(eq2_subs, Rn).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Nn,N_homo)])

    a3 = diff(eq3_subs, Nn).subs([(N,N_homo), (L,L_homo),             (Ln,L_homo)])
    b3 = diff(eq3_subs, Ln).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Ln,L_homo)])
    c3 = diff(eq3_subs, Rn).subs([(N,N_homo), (L,L_homo), (R,R_homo), (Ln,L_homo)])
    
    #And store these values into the 'B' array
    B = np.array([[a1, b1, c1], [a2, b2, c2], [a3, b3, c3]], dtype = 'float64')
    
    # The final matrix is H + qk*B, where qk is equal to -1 for an unidimensional arrangement of cells
    mat = H - B
    w, v = LA.eig(mat)
    
    return w.max() # The function returns the highest eigenvalue of the calculated matrix. 
# This value is the Maxium Lyapunov Exponent used to generate the final plot
```


```python
# Now we are going to generate the matrix that contains the MLE values depending on the kc and kt.
# This is the matrix 

n_points = 101
values = 10**(np.linspace(-3, 2, n_points))

# We generate matrix initially filled with zeros, were we are going to store the MLE values.
result = np.zeros((n_points,n_points))

# We initialize the parameters and variables in the equations

N, L, R, n, m, kc, kt, gamma, gamma_R, k_RS, beta_N,  beta_L, beta_R = symbols('N L R n m kc kt gamma gamma_R k_RS beta_N beta_L beta_R')#, real = True)

eq1 = beta_N - gamma*N - kt*N*L - kc*N*L
eq2 = beta_L*(1/(1+R**m)) - gamma*L - kt*N*L - kc*N*L
eq3 = beta_R*((N*L)**n)/(k_RS + (N*L)**n) - gamma_R*R

n_val = 1
m_val = 1
gamma_val = 1
gamma_R_val = 1
k_RS_val = 1
beta_N_val = 2
beta_L_val = 10
beta_R_val = 10

# We do a parameter scan for kc and kt with values logarithmically spaced
for i, kc_val in enumerate(values):
    print(kc_val)
    for j, kt_val in enumerate(values):
        eq1_subs = eq1.subs([(beta_N,beta_N_val),(gamma,gamma_val),(kt,kt_val),(kc,kc_val)])
        eq2_subs = eq2.subs([(beta_L,beta_L_val),(gamma,gamma_val),(kt,kt_val),(kc,kc_val), (m, m_val)])
        eq3_subs = eq3.subs([(beta_R,beta_R_val),(gamma_R,gamma_R_val),(k_RS, k_RS_val), (n, n_val)])

        def f(variables) :
            (N, L, R) = variables
            n = 1
            m = 1
            gamma = 1
            gamma_R = 1
            k_RS = 1
            beta_N = 2
            beta_L = 10
            beta_R = 10
            kc = kc_val
            kt = kt_val
        
            eq1 = beta_N - gamma*N - kt*N*L - kc*N*L
            eq2 = beta_L*(1/(1+R**m)) - gamma*L - kt*N*L - kc*N*L
            eq3 = beta_R*((N*L)**n)/(k_RS + (N*L)**n) - gamma_R*R

            return [eq1, eq2, eq3]
        
        # Here we calculate the values of N, L and R that correspond to the steady state
        steady_state = opt.fsolve(f, (0.5, 1.0, 1.0) )
        
        #With the following function we calculate the MLE value for a given steady state
        result[i,j] = calculate_MLE(steady_state, n_val = n_val, m_val = m_val, kc_val = kc_val, kt_val = kt_val,
                                                 gamma_val = gamma_val, gamma_R_val = gamma_R_val, k_RS_val = k_RS_val, 
                                                 beta_N_val = beta_N_val, beta_L_val = beta_L_val, beta_R_val = beta_R_val)
        
```


```python
# Here we define the position of the labels in the final plot
tick_values = [value for i,value in enumerate(values) if i%20 == 0]
tick_positions = np.arange(len(tick_values))*20 +0.5
```


```python
# And finally, we generate the final plot
plt.figure(figsize = (15,10))
sns.heatmap(np.flipud(result), cmap = 'seismic', center = 0.0)
sns.lineplot([0.001, 100], 101-np.array([0.001, 100]), color = '#FF00FF')
plt.xticks(tick_positions, tick_values)
plt.xlabel('kt')
plt.yticks(np.flip(tick_positions), tick_values)
plt.ylabel('kc')
```

9e3bf7934fed9beb9f1a050456e07de783664d09
ipynb
Jupyter Notebook
LSA/LSA.ipynb
hsancheziranzo/notochord-lateral-inhibition
2e5c5fe15e30ea6bacdc0282e1506b44b05415af
MIT
LSA/LSA.ipynb
hsancheziranzo/notochord-lateral-inhibition
2e5c5fe15e30ea6bacdc0282e1506b44b05415af
MIT
LSA/LSA.ipynb
hsancheziranzo/notochord-lateral-inhibition
2e5c5fe15e30ea6bacdc0282e1506b44b05415af
MIT
Qwen/Qwen-72B
1. YES
2. YES