<jupyter_start><jupyter_text>Lecture 2<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import sympy as sy
import tellurium as te
FONTSIZE  = 18<jupyter_output><empty_output><jupyter_text>Phase Portrait<jupyter_code>def plotPhase(xlim, ylim, dFunc, numPt=20):
    """
    Constructs a phase plot for a system of 2-d system of differential equations.

    Parameters
    ----------
    xlim: tuple
        lower and upper value of x
    ylim: tuple
        lower and upper value of y
    dFunc: Function
        args:
            tuple of current x, y
    numPt: int
        number of points in the grid
    """
    # Create the mesh
    xGrid = np.linspace(xlim[0], xlim[1], numPt)
    yGrid = np.linspace(ylim[0], ylim[1], numPt)
    xMesh, yMesh = np.meshgrid(xGrid, yGrid)
    xVals, yVals = np.zeros(xMesh.shape), np.zeros(yMesh.shape)
    # Populate mesh with values
    numX, numY = yMesh.shape
    for xIdx in range(numX):
        for yIdx in range(numY):
            x = xMesh[xIdx, yIdx]
            y = yMesh[xIdx, yIdx]
            yprime = dFunc([x, y])
            xVals[xIdx, yIdx] = yprime[0]
            yVals[xIdx, yIdx] = yprime[1]
    # Construct quiver plot
    fig, ax = plt.subplots(1, figsize=(10, 12))
    quiverPlt = ax.quiver(xMesh, yMesh, xVals, yVals, color='r')
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    
def dFunc(Y):
    y1, y2 = Y
    return [y2, -np.sin(y1)]
numPt = 50
xlim = [-2, 8]
ylim = [-4, 4]

plotPhase(xlim, ylim, dFunc, numPt=numPt)<jupyter_output><empty_output><jupyter_text>1-dimensional System $\frac{dS_1}{dt} = -(k_f +k_r)S_1 + k_r S_T$<jupyter_code>C_k_f = 1
C_k_r = 0.8
C_S_T = 10
model = """
S1 -> S2; k_f*S1
S2 -> S1; k_r*S2

k_f = 0
k_r = 0
S1 = 0
S2 = 0
"""
rr = te.loada(model)
rr.k_f = C_k_f
rr.k_r = C_k_r
rr.S1 = C_S_T
data = rr.simulate()
fig, ax = plt.subplots(figsize=(5, 5))
_ = ax.plot(data[:,0], data[:, 1:])
ax.set_ylabel("concentration", fontsize=FONTSIZE)
ax.set_xlabel("time (sec)", fontsize=FONTSIZE)
_ = ax.legend(data.colnames[1:])
values = np.array([0.1*v for v in range(0, 10*C_S_T+1)])
def ds(s):
     return np.array(-(C_k_f + C_k_r)*s + C_k_r*C_S_T)
fig, ax = plt.subplots(1, figsize=(8, 6))
ax.plot(values, ds(values))
ax.set_xlabel(r'$S_1$', fontsize=FONTSIZE)
ax.set_ylabel(r'$\frac{dS_1}{dt}$', fontsize=FONTSIZE, rotation=0)
_ = ax.plot([0, 10], [0, 0], linestyle="--")<jupyter_output><empty_output><jupyter_text>$\frac{dS_1}{dt} = -(k_f +k_r)S_1(S_T - S_1) + k_r S_T = $<jupyter_code>sy.init_printing()
dS1, S1, k_f, k_r, S_T, k_2, k_3, k_0, S2 = sy.symbols("dS1 S1 k_f k_r S_T k_2 k_3 k_0 S2")
f_dS1 = -(k_f + k_r)*S1*(S_T - S1) + (k_r)*S_T
sy.simplify(f_dS1)
ans = sy.solve(f_dS1, S1)
ans
ans[0].evalf(subs={S_T: C_S_T, k_f: C_k_f, k_r: C_k_r})<jupyter_output><empty_output><jupyter_text>2-d Linear System- $\emptyset \rightarrow S_1$; $k_0$- $S_1 \leftrightarrow S_2$; $-k_f S_1 + k_r S_2$- $S_1 \rightarrow 2 S_1$; $k_2 S_1$- $S_2 \rightarrow \emptyset$; $k_3 S_2$$\frac{S_1}{dt} = k_0 -k_f S_1 + k_r S_2 + k_2 S_1 = k_0 + S_1 (k_2 - k_f) + k_r S_2$$\frac{S_2}{dt} = k_f S_1 - k_r S_2 - k_3 S_2 = k_f S_1 - S_2 (k_r + k_3)$ The model turns unstable when k2 transitions from .5 to .6.<jupyter_code>model = """
J0:  -> S1; k_0
J1: S1 -> S2; k_f*S1
J2: S2 -> S1; k_r*S2
J3: S1 -> 2 S1; k_2*S1
J4: S2 -> ; k_3*S2

k_0 = 1
k_2 = 0.5
k_3 = 1
k_f = 1.0
k_r = 0.8
S1 = 180
S2 = 100

//at (time > 10): S1 = 200
"""
rr = te.loada(model)
#rr.k_f = C_k_f
#rr.k_r = C_k_r
#rr.S_T = C_S_T
data = rr.simulate(0, 100)
rr.plot(data)
dS1 = k_0 + S1*(k_2 - k_f) + k_r*S2
dS2 = k_f*S1 - S2*(k_3 + k_r)
C_k_f = 1
C_k_r = 0.8
C_k_0 = 1
C_k_2 = .55
C_k_3 = 1
def dFunc(Svec):
    s1, s2 = Svec
    ds1 = dS1.subs({k_0: C_k_0, S1: s1, S2: s2, k_2: C_k_2, k_f: C_k_f, k_r: C_k_r})
    ds2 = dS2.subs({k_0: C_k_0, S1: s1, S2: s2, k_2: C_k_2, k_f: C_k_f, k_r: C_k_r, k_3: C_k_3})
    return (ds1, ds2)
subDct = {k_0: C_k_0, k_2: C_k_2, k_f: C_k_f, k_r: C_k_r, k_3: C_k_3}
def findFixedPoints(expressions, subDct, vars):
    newExpressions = []
    for expression in expressions:
        newExpressions.append(expression.subs(subDct))
    return sy.solve(newExpressions, vars)

findFixedPoints([dS1, dS2], subDct, [S1, S2])
plotPhase([0, 50], [0, 50], dFunc, numPt=20)
xx = k_0 + S1
xx.subs({k_0: 1})
dS1.subs({k_0:1, k_2:1, k_f:1, k_r: 1})
dS1.find(k_0)
vMat = sy.Matrix([dS1, dS2])
vMat
jacobian = vMat.jacobian((S1, S2))
jacobian
det = sy.simplify(jacobian.det())
det<jupyter_output><empty_output><jupyter_text>Note that if $k_2 = 0$, then the determinant is always positive.<jupyter_code>fixedPt = sy.solve([dS1, dS2], [S1, S2])
fixedPt<jupyter_output><empty_output><jupyter_text>Lac Operon<jupyter_code>dI_int, a, I_ext, I_int, P, b, c, e, dP, I_int_str, P_str = sy.symbols("dI_int a I_ext I_int P b c e dP I_int_str P_str")
dP = c*I_int - e*P
dI_int = a*(I_ext - I_int)*P - b*I_int
fixedPoints = sy.solve([dP, dI_int], [I_int, P])
fixedPoints[1]
subDct = {a:1, b:0.1, c:1, e:0.5, I_ext:3}
C_I_int_str, C_P_str = [s.evalf(subs=subDct) for s in fixedPoints[1]]
subDct[I_int_str] = C_I_int_str
subDct[P_str] = C_P_str
# Jacobian
jacobian = sy.Matrix([dI_int, dP])
jacobian = jacobian.jacobian([I_int, P])
linearization = jacobian.subs({P: P_str, I_int:I_int_str})
linearization.subs(subDct)
2.95*4<jupyter_output><empty_output><jupyter_text>Nonlinear System<jupyter_code>values = np.array([0.1*v for v in range(0, 10*C_S_T+1)])
def ds(s):
     return np.array(-(C_k_f + C_k_r)*s*(C_S_T-s) + C_k_r*C_S_T)
fig, ax = plt.subplots(1, figsize=(8, 6))
ax.plot(values, ds(values))
ax.set_xlabel(r'$S_1$', fontsize=FONTSIZE)
ax.set_ylabel(r'$\frac{dS_1}{dt}$', fontsize=FONTSIZE, rotation=0)
_ = ax.plot([0, 10], [0, 0], linestyle="--")
dx3, dx1, x1, ll, k1, y1, y2, x2, dx2, k2 = sy.symbols("dx3, dx1, x1, ll, k1, y1, y2, x2, dx2, k2 ")
dx1 = ll - k1*y1*(x1/(1 + x1))
dx2 = k1*y1*(x1/(1 + x1)) - k2*y2*(x2/1+x2)
fixedPoints = sy.solve([dx1, dx2], [x1, x2])
fixedPoints
subDct = {x2: fixedPoints[0][1]}
dx3 = k2*y2*(x2/1+x2)
dx3.evalf(subs=subDct)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1768562530732667, 'fraction_numerical': 0.05622029175544992, 'mean_word_length': 3.920967741935484, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '643130', 'license': 'MIT', 'n_tokens_mistral': 3220, 'n_tokens_neox': 3118, 'n_words': 647, 'path': 'archived_lectures/Winter 2021/lecture_2/lecture_2.ipynb', 'repo_name': 'BioModelTools/topics-course', 'size': '6101'}