In this notebook we will look at how we can explore a simple kinetic system using only a handful of lines of Python code and the ChemPy library. First we import the functions & classes that we will be using:
from collections import defaultdict
from ipywidgets import interact, FloatSlider
from chempy import ReactionSystem, Substance
from chempy.kinetics.ode import get_odesys
from chempy.kinetics._native import get_native
from chempy.units import SI_base_registry, default_units as u
import matplotlib.pyplot as plt
%matplotlib inline
Next we define our what reactions are taking place:
rsys = ReactionSystem.from_string("""
A -> B; 'k1'
B + C -> P; 'k2'
""", substance_factory=Substance)
rsys
From this reactionsystem we generate a system of ordinary differential equations:
odesys, extra = get_odesys(rsys, include_params=False, unit_registry=SI_base_registry,
output_conc_unit=u.micromolar, output_time_unit=u.minute)
Not that it is strictly needed for this simple system, but since we are using pyodesys, we can generate C++ code which is compiled to a native extension module:
native = get_native(rsys, odesys, 'cvode')
We will use ipywidgets for our interactive controls:
def integrate_and_plot(
tend_minutes=FloatSlider(5, min=1, max=60, step=1),
lgA0molar=FloatSlider(-6, min=-7, max=-5, step=.1),
lgC0molar=FloatSlider(-6, min=-7, max=-5, step=.1),
k1_per_min=FloatSlider(5.8, min=1.0, max=20.0, step=0.1),
k2_per_M_per_s=FloatSlider(4000, min=500, max=16000, step=500)
):
result = native.integrate(tend_minutes*u.minute, defaultdict(lambda: 0*u.M, {'A': 10**lgA0molar*u.M, 'C': 10**lgC0molar*u.M}),
{'k1': k1_per_min/u.minute, 'k2': k2_per_M_per_s/u.M/u.s}, integrator='cvode')
result.plot(title_info=2)
Try exploring the system by adjusting the parameter sliders below:
interact(integrate_and_plot)