Solving Ordinary Differential Equations (ODEs)

import numpy as np

from scipy.integrate import odeint

import matplotlib.pyplot as plt

%matplotlib inline

# function that returns dy/dt

def model(y,t,k):

    dydt = -k * y

    return dydt

# initial condition

y0 = 5

# time points

t = np.linspace(0, 20)

# solve ODEs

y1 = odeint(model,y0,t,args=(0.1,))

y2 = odeint(model,y0,t,args=(0.2,))

y3 = odeint(model,y0,t,args=(0.5,))

# plot results

plt.figure()

plt.plot(t,y1,'r-',linewidth=2,label='k=0.1')

plt.plot(t,y2,'b--',linewidth=2,label='k=0.2')

plt.plot(t,y3,'g:',linewidth=2,label='k=0.5')

plt.xlabel('time')

plt.ylabel('y(t)')

plt.legend()

"""

Hodgkin-Huxley model of excitable barnacle muscle fiber

reviewed in Rinzel (1990) Bulletin of Mathematical Biology 52 pp. 5-23.

"""

import numpy as np

from scipy.integrate import odeint

import matplotlib.pyplot as plt

from math import exp, expm1

def get_iStim(t):

    if 20 < t <= 21:

        return -6.65

    elif 60 < t <= 61:

        return -6.86

    else:

        return 0

# Model

def hh_rhs(y, t, p):

    # Constants

    E_N = p['E_N']  # Reversal potential of Na

    E_K = p['E_K']  # Reversal potential of K

    E_LEAK = p['E_LEAK']  # Reversal potential of leaky channels

    G_N_BAR = p['G_N_BAR']  # Max. Na channel conductance

    G_K_BAR = p['G_K_BAR'] # Max. K channel conductance

    G_LEAK = p['G_LEAK']  # Max. leak channel conductance

    C_M = p['C_M']  # membrane capacitance

    # Equations

    v, m, h, n = y

    mAlfaV = -0.10 * (v + 35)

    mAlfa = mAlfaV / expm1(mAlfaV)

    mBeta = 4.0 * exp(-(v + 60) / 18.0)

    dm = -(mAlfa + mBeta) * m + mAlfa

    hAlfa = 0.07 * exp(-(v+60)/20)

    hBeta = 1 / (exp(-(v+30)/10) + 1)

    dh  = -(hAlfa + hBeta) * h + hAlfa

    iNa = G_N_BAR * (v - E_N) * (m**3) * h

    nAlfaV = -0.1 * (v+50)

    nAlfa = 0.1 * nAlfaV / expm1(nAlfaV)

    nBeta = 0.125 * exp( -(v+60) / 80)

    dn  = -(nAlfa + nBeta) * n + nAlfa

    iK = G_K_BAR * (v - E_K) * (n**4)

    iLeak = G_LEAK * (v - E_LEAK)

    iSt = get_iStim(t)

    dv = -(iNa + iK + iLeak + iSt) / C_M

    return [dv, dm, dh, dn]

# Initial conditions

y0 = v, m, h, n = -59.8977, 0.0536, 0.5925, 0.3192

# Parameters

p = {'E_N': 55,     # Reversal potential of Na

     'E_K': -72,    # Reversal potential of K

     'E_LEAK': -49, # Reversal potential of leaky channels

     'G_N_BAR': 120,# Max. Na channel conductance

     'G_K_BAR': 36, # Max. K channel conductance

     'G_LEAK': 0.3, # Max. leak channel conductance

     'C_M': 1.0}    # membrane capacitance

# time span

tStart, tEnd = 0, 100

N_POINTS = 100+1

ts = np.linspace(tStart, tEnd, N_POINTS)

# Solve the ODEs

sol = odeint(hh_rhs, y0, ts, args=(p,), tcrit=[20, 21, 60, 61])

# Plotting

fig, axs = plt.subplots(nrows=2, figsize=(10, 15))

axs[0].plot(ts, sol[:, 0], 'k-', linewidth=2)

axs[0].set_xlabel("Time (ms)")

axs[0].set_ylabel("Membrane voltage (mV)")

axs[1].plot(ts, sol[:, 1:], linewidth=2, label=["m", "h", "n"])

axs[1].set_xlabel("Time (ms)")

axs[1].set_ylabel("Gating variables")

fig

conda install -c alubbock pysb kappa --yes

from pysb import *

from pysb.simulator import ScipyOdeSimulator

from pylab import linspace, plot, xlabel, ylabel, gcf

import matplotlib.pyplot as plt

# A simple model with a reversible binding rule

Model()

# Declare the monomers

Monomer('L', ['s'])

Monomer('R', ['s'])

# Declare the parameters

Parameter('L_0', 100)

Parameter('R_0', 200)

Parameter('kf', 1e-3)

Parameter('kr', 1e-3)

# Declare the initial conditions

Initial(L(s=None), L_0)

Initial(R(s=None), R_0)

# Declare the binding rule

Rule('L_binds_R', L(s=None) + R(s=None) | L(s=1) % R(s=1), kf, kr)

# Observe the complex

Observable('LR', L(s=1) % R(s=1))

time = linspace(0, 40, 100)

print("Simulating...")

sim_result = ScipyOdeSimulator(model, time).run()

# Plot the trajectory of LR

plt.figure()

plot(time, sim_result.observables['LR'])

xlabel('Time (seconds)')

ylabel('Amount of LR')

gcf()

from pysb import Model, Monomer, Parameter, Initial, Rule, Observable

from pysb.macros import *

# Create an instance of empty model

Model()

# Create a monomer named 'enzyme' that only has one attribute, called 'binding1'.

Monomer('enzyme', ['binding1'])

# Create a monomer named 'protein' that has two attributes, called 'binding1'(bound to enzyme or not) and state (substrate or product).

Monomer('protein', ['binding', 'state'], {'state': ['sub','pro']})

# Define a forward, reverse and catalysis parameters (kf, kr, kc)

Parameter('kf',0.003)

Parameter('kr',0.001)

Parameter('kc',0.002)

# Define a binding rule where the unbound enzyme binds the protein in the substrate state. | = reversible reaction, % = bounded state

Rule('binding', enzyme(binding1=None) + protein(state='sub', binding=None)

     | enzyme(binding1=1) % protein(state='sub', binding=1), kf, kr)

# Create a catalysis rule where the complex dissociates and the protein change state to product (the enzyme is left unaltered). >>  = irreversible reaction

Rule('dissociation', enzyme(binding1=1) % protein(state='sub', binding=1) 

     >> enzyme(binding1=None) + protein(state='pro', binding=None), kc)

# Define the initial conditions for the two objects using parameters.

Parameter('enzyme_0', 100)

Parameter('protein_0', 50)

Initial(enzyme(binding1=None), enzyme_0 )

Initial(protein(binding=None, state='sub') , protein_0)

# Define the observables: total enzyme, product and the complex

Observable('e_total', enzyme())

Observable('e_free', enzyme(binding1=None))

Observable('substrate', protein(binding=None, state='sub'))

Observable('product', protein(binding=None, state='pro'))

Observable('complex', enzyme(binding1=1) % 

           protein(binding=1, state='sub'))

# run a simulation

from pysb.integrate import odesolve

import matplotlib.pyplot as plt

import numpy as np

tspan = np.linspace(0, 100, 100)

y = odesolve(model, tspan)

plt.figure()

plt.plot(tspan, y['substrate'], label="substrate")

plt.plot(tspan, y['e_total'], label="e_total")

plt.plot(tspan, y['e_free'], label="e_free")

plt.plot(tspan, y['product'], label="product")

plt.plot(tspan, y['complex'], label="complex")

plt.xlabel('time')

plt.ylabel('population')

plt.ylim(0,110)

plt.legend(loc='best')

plt.gcf()