chemistrykit.kinetics#
chemistrykit.kinetics: reaction kinetics.
Elementary integrated rate laws and half-lives, the Arrhenius equation and
activation-energy fitting, Michaelis-Menten enzyme kinetics with
Lineweaver-Burk linearization and inhibition, a general
stoichiometry-matrix reaction-network engine (parallel, consecutive,
reversible, and steady-state-approximation chain reactions) integrated
via chemistrykit.integrators, the Brusselator oscillating
reaction network and the Oregonator model of the Belousov-Zhabotinsky
reaction, collision/diffusion/transition-state theories of the rate
constant, and Gillespie’s exact stochastic simulation algorithm.
- class chemistrykit.kinetics.ArrheniusFit(Ea, A, r_squared, R_gas=8.31446261815324)[source]#
Bases:
objectResult of fitting rate-constant-vs-temperature data to the Arrhenius equation.
- class chemistrykit.kinetics.Brusselator(X0=1.0, Y0=1.0, A=1.0, B=3.0)[source]#
Bases:
ReactionNetworkThe Brusselator oscillating reaction network.
Derived from four elementary steps (with B, D, E, and the source of X held as constant reservoirs; Prigogine & Lefever 1968):
\[\begin{split}A &\to X \\ 2X + Y &\to 3X \\ B + X &\to Y + D \\ X &\to E\end{split}\]which, in units where all four rate constants are 1, gives the net kinetics
\[\frac{dX}{dt} = A - (B+1)X + X^2 Y, \qquad \frac{dY}{dt} = BX - X^2 Y\]The single fixed point \((X^*, Y^*) = (A, B/A)\) (see
fixed_point()) undergoes a Hopf bifurcation at \(B = 1 + A^2\) (seeis_above_hopf_threshold()): below threshold the fixed point is a stable focus and concentrations relax to it; above threshold it is an unstable focus surrounded by a stable limit cycle, so concentrations settle into sustained oscillation instead.- Parameters:
- fixed_point()[source]#
Return the network’s single steady state
(X*, Y*) = (A, B/A).- Return type:
- Returns:
ndarray, shape (2,)
Examples
The fixed point is exactly a zero of the vector field:
>>> import numpy as np >>> system = Brusselator(A=1.0, B=3.0) >>> fp = system.fixed_point() >>> np.allclose(system.rhs(fp, 0.0), 0.0) True
- is_above_hopf_threshold()[source]#
Whether
B > 1 + A^2, the Hopf-bifurcation onset of sustained oscillation.- Return type:
- Returns:
bool
Examples
>>> Brusselator(A=1.0, B=3.0).is_above_hopf_threshold() True >>> Brusselator(A=1.0, B=1.5).is_above_hopf_threshold() False
- class chemistrykit.kinetics.EyringFit(dH, dS, r_squared)[source]#
Bases:
objectResult of fitting rate-constant-vs-temperature data to the Eyring equation.
- class chemistrykit.kinetics.FirstOrder(k, C0)[source]#
Bases:
RateLawFirst-order kinetics:
d[A]/dt = -k[A].Integrated law: \([A](t) = [A]_0 e^{-kt}\). The half-life \(t_{1/2} = \ln(2)/k\) is independent of the initial concentration – the defining signature of first-order kinetics (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 20, Table 20.3), used e.g. for radioactive decay and many unimolecular reactions.
- Parameters:
Examples
>>> import numpy as np >>> law = FirstOrder(k=np.log(2), C0=1.0) >>> round(float(law.concentration(1.0)), 6) 0.5 >>> round(law.half_life(), 6) 1.0
- class chemistrykit.kinetics.KineticsResult(t, y, species=<factory>, method='', extra=<factory>)[source]#
Bases:
objectContainer for the output of a
ReactionNetwork.integrate()call.Mirrors physicskit’s
SimulationResult(a stable, dataclass return type consumed by visualizers), specialized to a named-species concentration trajectory instead of a phase-space one.- concentration(name)[source]#
Return the concentration trajectory of a single named species.
- Parameters:
name (
str) – A species name fromself.species.- Return type:
- Returns:
ndarray, shape (n_steps + 1,)
Examples
>>> import numpy as np >>> result = KineticsResult(t=np.array([0.0, 1.0]), y=np.array([[1.0, 0.0], [0.5, 0.5]]), species=("A", "B")) >>> result.concentration("B") array([0. , 0.5])
- class chemistrykit.kinetics.MichaelisMentenFit(Vmax, Km, r_squared)[source]#
Bases:
objectResult of a Lineweaver-Burk (double-reciprocal) linear fit.
- class chemistrykit.kinetics.MichaelisMentenProgress(S0, Vmax, Km)[source]#
Bases:
ReactionNetworkSubstrate-depletion progress curve under Michaelis-Menten kinetics.
Integrates
d[S]/dt = -Vmax*[S]/(Km+[S])forward in time viachemistrykit.integrators– the “progress curve” an enzyme assay actually measures, as distinct from the initial-ratemichaelis_menten_rate()used for a Lineweaver-Burk fit (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 20.6). In the two asymptotic limits this reduces to the closed-form rate laws inchemistrykit.kinetics.systems.rate_laws: zeroth order in S whenS0 >> Km(ZeroOrder) and first order whenS0 << Km(FirstOrder).- Parameters:
- class chemistrykit.kinetics.Oregonator(x0=0.5, y0=0.1, z0=0.1, epsilon=0.04, epsilon_prime=0.0004, q=0.0008, f=1.0)[source]#
Bases:
ReactionNetworkThe Field-Noyes Oregonator model of the Belousov-Zhabotinsky reaction.
Field and Noyes (J. Chem. Phys. 60, 1877 (1974)) reduced the FKN mechanism of the cerium-catalyzed bromate/malonic-acid (BZ) reaction to five steps among three intermediates – X = HBrO2, Y = Br-, and Z = Ce(IV) – with the bromate (A) and organic substrate (B) pools held constant:
\[\begin{split}A + Y &\to X + P, \qquad X + Y \to 2P, \qquad A + X \to 2X + 2Z, \\ 2X &\to A + P, \qquad B + Z \to f\,Y\end{split}\]In the standard dimensionless scaling (Tyson, in R. J. Field and M. Burger (eds.), Oscillations and Traveling Waves in Chemical Systems, Wiley 1985, Ch. 3) the kinetics read
\[\epsilon \frac{dx}{dt} = qy - xy + x(1-x), \quad \epsilon' \frac{dy}{dt} = -qy - xy + fz, \quad \frac{dz}{dt} = x - z\]The system is stiff (\(\epsilon' \ll \epsilon \ll 1\)), so integrate with
method="dopri5". For the default parameters it relaxes onto a stable limit cycle (sustained BZ oscillation); for large stoichiometric factor f (e.g.f=3) its steady state is stable instead.- Parameters:
x0 (
float) – Initial scaled concentrations.y0 (
float) – Initial scaled concentrations.z0 (
float) – Initial scaled concentrations.epsilon (
float) – Scaled Oregonator parameters (all positive).epsilon_prime (
float) – Scaled Oregonator parameters (all positive).q (
float) – Scaled Oregonator parameters (all positive).f (
float) – Scaled Oregonator parameters (all positive).
- fixed_point()[source]#
Return the positive steady state
(x*, y*, z*).Setting all three derivatives to zero gives \(z^* = x^*\), \(y^* = f x^*/(q + x^*)\), and \(x^*\) as the positive root of \(x^2 - (1 - f - q)x - q(1 + f) = 0\).
- Return type:
- Returns:
ndarray, shape (3,)
Examples
>>> import numpy as np >>> system = Oregonator(f=1.0) >>> np.allclose(system.rhs(system.fixed_point(), 0.0), 0.0, atol=1e-9) True
- class chemistrykit.kinetics.RateLaw[source]#
Bases:
ABCBase class for a single elementary reaction’s integrated rate law.
Subclasses implement the closed-form solution of
dC/dt = -rate(C)(or, for zero order,+ratewith a sign convention appropriate to that order) directly – there is nothing to numerically integrate for these textbook cases, unlikeReactionNetwork.
- class chemistrykit.kinetics.ReactionNetwork(state0)[source]#
Bases:
ABCCommon base for a system of species evolving as
dC/dt = f(C, t).Concrete subclasses must set, in
__init__:self.species: tuple of str – ordered species names.self._rhs_njit: an@njitdispatcher with signature(state, t, params) -> dstate(thechemistrykit.integrators.RHSFuncconvention).self.params: ndarray – the parameter vector read byself._rhs_njit(ornp.empty(0)if every rate constant is instead baked into the njit closure itself, aschemistrykit.kinetics.systems.networks’s mass-action engine does with its stoichiometry/rate-constant arrays).
and call
super().__init__(state0).- integrate(t_span, dt=None, method='rk4', **kwargs)[source]#
Integrate the reaction network forward in time.
- Parameters:
dt (float, optional) – Fixed step size, required for
method="rk4". Formethod="dopri5"it is instead used as the initial step size attempt (defaulting to(t1 - t0) / 1000), since the adaptive integrator adjusts it automatically.method (
str) – Integrator to use."dopri5"(adaptive Dormand-Prince) is recommended for stiff-ish networks with widely separated rate constants (e.g. a steady-state-approximation chain).**kwargs – Forwarded to
chemistrykit.integrators.dopri5_integrate()(e.g.rtol,atol) whenmethod="dopri5".
- Return type:
- Returns:
KineticsResult
- reset(state0=None, t0=0.0)[source]#
Reset the system’s state and clock.
- Parameters:
state0 (array-like, optional) – New state; if omitted, the current state is kept.
t0 (
float) – New time.
- Returns:
ndarray – The (possibly updated) current state.
- class chemistrykit.kinetics.SecondOrder(k, C0)[source]#
Bases:
RateLawSecond-order kinetics (single reactant):
d[A]/dt = -k[A]^2.Integrated law: \(1/[A](t) = 1/[A]_0 + kt\), i.e. \([A](t) = [A]_0 / (1 + k[A]_0 t)\). Unlike first order, the half-life \(t_{1/2} = 1/(k[A]_0)\) does depend on the initial concentration (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 20, Table 20.3).
- Parameters:
Examples
>>> law = SecondOrder(k=1.0, C0=1.0) >>> round(float(law.concentration(1.0)), 6) 0.5 >>> round(law.half_life(), 6) 1.0
- class chemistrykit.kinetics.StochasticTrajectory(t, counts, species=<factory>)[source]#
Bases:
objectOne exact sample path of a stochastic reaction network.
- counts: ndarray#
Molecule counts immediately after each event (row 0 is the initial state).
- Type:
ndarray of int, shape (n_events + 1, n_species)
- class chemistrykit.kinetics.StoichiometricNetwork(species, stoich_matrix, rate_constants, reactant_orders, state0)[source]#
Bases:
ReactionNetworkA network of elementary mass-action reactions, integrated numerically.
- Parameters:
stoich_matrix (array-like, shape (n_species, n_reactions)) – Net stoichiometric coefficient of each species in each reaction (negative for a net-consumed species, positive for net-produced).
rate_constants (array-like, shape (n_reactions,)) – Rate constant of each reaction.
reactant_orders (array-like, shape (n_species, n_reactions)) – Reactant-side order of each species in each reaction’s mass-action rate law (0 for a species that is not a reactant in that reaction; for an elementary reaction this equals the reactant’s stoichiometric coefficient, but is given independently so non-elementary/pseudo-order rate laws can also be modeled).
state0 (array-like, shape (n_species,)) – Initial concentrations, in species order.
- classmethod consecutive(k1, k2, A0=1.0)[source]#
Build the consecutive (chain) mechanism
A -> B(k1),B -> C(k2).- Parameters:
- Return type:
- Returns:
StoichiometricNetwork
Examples
>>> net = StoichiometricNetwork.consecutive(k1=1.0, k2=0.3, A0=1.0) >>> result = net.integrate((0.0, 20.0), dt=1e-3, method="rk4") >>> round(float(result.concentration("A")[-1] + result.concentration("B")[-1] + result.concentration("C")[-1]), 6) 1.0
- classmethod parallel(k1, k2, A0=1.0)[source]#
Build the parallel (competing) mechanism
A -> B(k1),A -> C(k2).- Parameters:
- Return type:
- Returns:
StoichiometricNetwork
Examples
The product ratio [B]/[C] is exactly k1/k2 at every time, the classic parallel-reaction result (Espenson, Chemical Kinetics and Reaction Mechanisms, 2nd ed., Ch. 4):
>>> import numpy as np >>> net = StoichiometricNetwork.parallel(k1=2.0, k2=1.0, A0=1.0) >>> result = net.integrate((0.0, 5.0), dt=1e-3, method="rk4") >>> round(float(result.concentration("B")[-1] / result.concentration("C")[-1]), 4) 2.0
- classmethod reversible(kf, kr, A0=1.0, B0=0.0)[source]#
Build the reversible mechanism
A <-> B(forward kf, reverse kr).- Parameters:
- Return type:
- Returns:
StoichiometricNetwork
Examples
>>> net = StoichiometricNetwork.reversible(kf=2.0, kr=1.0, A0=1.0) >>> result = net.integrate((0.0, 20.0), dt=1e-3, method="rk4") >>> round(float(result.concentration("A")[-1]), 3) 0.333
- class chemistrykit.kinetics.ZeroOrder(k, C0)[source]#
Bases:
RateLawZero-order kinetics:
d[A]/dt = -k.Integrated law: \([A](t) = [A]_0 - kt\), valid only while \([A](t) \geq 0\) – physically, the reaction stops (or crosses over to a different rate law) once the reactant is exhausted, at \(t = [A]_0 / k\).
- Parameters:
Examples
>>> law = ZeroOrder(k=0.1, C0=1.0) >>> round(float(law.concentration(5.0)), 6) 0.5 >>> round(law.half_life(), 6) 5.0
- rate(t=None)[source]#
Return the (constant) rate
k.- Parameters:
t (float or array-like of float, optional) – Unused (the rate does not depend on time or concentration for a zero-order reaction); accepted for interface consistency with
FirstOrder/SecondOrder.- Returns:
float
- chemistrykit.kinetics.arrhenius_rate_constant(A, Ea, T, R_gas=8.31446261815324)[source]#
Evaluate the Arrhenius equation \(k = A e^{-E_a/(RT)}\).
- Parameters:
- Returns:
float or ndarray – Rate constant(s) k, in the same units as A.
Examples
>>> round(float(arrhenius_rate_constant(A=1e13, Ea=50e3, T=300.0)), 2) 19696.84
- chemistrykit.kinetics.collision_theory_rate_constant(T, sigma, reduced_mass, Ea=0.0, steric_factor=1.0)[source]#
Hard-sphere collision-theory rate constant of a bimolecular reaction.
Every collision between reactant molecules with relative kinetic energy along the line of centres above \(E_a\) reacts (times an empirical steric factor \(P\)):
\[k = P\,\sigma \sqrt{\frac{8 k_B T}{\pi \mu}}\; N_A\, e^{-E_a/RT}\]where \(\sigma = \pi d^2\) is the collision cross-section and \(\sqrt{8k_BT/\pi\mu}\) the mean relative speed (Atkins & de Paula, Physical Chemistry, 11th ed., Topic 18A).
- Parameters:
T (float or array-like of float) – Absolute temperature(s), in K.
sigma (
float) – Collision cross-section, in m^2.reduced_mass (
float) – Reduced mass \(\mu = m_A m_B/(m_A + m_B)\) of the colliding pair, in kg.Ea (
float) – Activation (threshold) energy, in J/mol.steric_factor (
float) – Steric factor \(P\).
- Returns:
float or ndarray – Rate constant, in m^3 mol^-1 s^-1.
Examples
With no barrier, \(k \propto \sqrt{T}\) exactly:
>>> import numpy as np >>> k = collision_theory_rate_constant(np.array([300.0, 1200.0]), sigma=4e-19, reduced_mass=1e-26) >>> round(float(k[1] / k[0]), 12) 2.0
- chemistrykit.kinetics.competitive_inhibition_rate(S, I, Vmax, Km, Ki)[source]#
Competitive-inhibition Michaelis-Menten rate law.
The inhibitor competes with substrate for the active site, which is equivalent to inflating the apparent \(K_m\) by a factor \(\alpha = 1 + [I]/K_i\) while leaving \(V_{max}\) unchanged (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 20.6):
\[v = \frac{V_{max}[S]}{K_m (1 + [I]/K_i) + [S]}\]- Parameters:
- Returns:
float or ndarray
Examples
With no inhibitor (
I=0) this reduces exactly to the uninhibited rate law:>>> float(competitive_inhibition_rate(S=2.0, I=0.0, Vmax=10.0, Km=2.0, Ki=1.0)) 5.0
- chemistrykit.kinetics.consecutive_analytic(A0, k1, k2, t)[source]#
Closed-form
A -> B -> Cconcentrations (Bateman’s equations).See Espenson, Chemical Kinetics and Reaction Mechanisms, 2nd ed., Ch. 4.2, or the original: H. Bateman, Proc. Cambridge Philos. Soc. 15, 423 (1910) (radioactive decay chains have the identical mathematical structure). The
k1 != k2formula below has a removable singularity atk1 == k2; the degenerate limit \([B](t) = A_0 k_1 t\,e^{-k_1 t}\) is used automatically in that case (obtained by L’Hopital’s rule / directly integrating \(d[B]/dt = k_1[A] - k_1[B]\) with equal rate constants).- Parameters:
- Returns:
A, B, C (float or ndarray)
Examples
Mass balance holds exactly at every time:
>>> import numpy as np >>> A, B, C = consecutive_analytic(A0=1.0, k1=1.0, k2=0.3, t=np.array([0.0, 1.0, 5.0, 20.0])) >>> np.allclose(A + B + C, 1.0) True
- chemistrykit.kinetics.diffusion_limited_rate_constant(T, viscosity)[source]#
Diffusion-limited rate constant from the solvent viscosity alone.
Combining
smoluchowski_rate_constant()with the Stokes-Einstein relation \(D = k_B T/(6\pi\eta r)\) for two reactants of equal radius \(r\) (so \(D = D_A + D_B = 2k_BT/(6\pi\eta r)\) and \(R^* = 2r\)) makes the radius cancel:\[k_D = \frac{8RT}{3\eta}\]- Parameters:
- Returns:
float or ndarray – Rate constant, in m^3 mol^-1 s^-1.
Examples
Water at 298 K (\(\eta \approx 8.9\times 10^{-4}\) Pa s) gives the familiar \(\sim 7\times 10^{9}\) L mol^-1 s^-1:
>>> k = diffusion_limited_rate_constant(298.15, 8.9e-4) >>> round(float(k) * 1000 / 1e9, 2) 7.43
- chemistrykit.kinetics.eyring_rate_constant(T, dH, dS, kappa=1.0)[source]#
Eyring (transition-state theory) rate constant.
\[k = \kappa \frac{k_B T}{h}\, e^{\Delta S^{\ddagger}/R}\, e^{-\Delta H^{\ddagger}/RT}\]i.e. \(k = \kappa (k_BT/h) e^{-\Delta G^{\ddagger}/RT}\) with \(\Delta G^{\ddagger} = \Delta H^{\ddagger} - T\Delta S^{\ddagger}\) (H. Eyring, J. Chem. Phys. 3, 107 (1935); Atkins & de Paula, Physical Chemistry, 11th ed., Topic 18C). Written for a unimolecular step (or with the standard-state concentration factor absorbed into \(\Delta S^{\ddagger}\)), so k is in s^-1.
- Parameters:
- Returns:
float or ndarray – Rate constant(s), in s^-1.
Examples
With zero activation enthalpy and entropy, k is the universal frequency \(k_BT/h\) (about 6.25e12 s^-1 at 300 K):
>>> round(float(eyring_rate_constant(300.0, dH=0.0, dS=0.0)) / 1e12, 3) 6.251
- chemistrykit.kinetics.fit_arrhenius(T, k, R_gas=8.31446261815324)[source]#
Fit rate constant vs. temperature data to the Arrhenius equation.
Linearizes \(\ln k = \ln A - E_a/R \cdot (1/T)\) and fits by ordinary least squares – the standard “Arrhenius plot” method (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 20.3): plotting \(\ln k\) against \(1/T\) gives a straight line of slope \(-E_a/R\) and intercept \(\ln A\).
- Parameters:
- Return type:
- Returns:
ArrheniusFit
Examples
Generate exact data from a known (Ea, A) and recover them:
>>> import numpy as np >>> T = np.array([280.0, 300.0, 320.0, 340.0, 360.0]) >>> k = arrhenius_rate_constant(A=5e12, Ea=60e3, T=T) >>> fit = fit_arrhenius(T, k) >>> round(fit.Ea, 2) 60000.0 >>> round(fit.r_squared, 6) 1.0
- chemistrykit.kinetics.fit_eyring(T, k)[source]#
Fit rate constant vs. temperature data to the Eyring equation.
Linearizes \(\ln(k/T) = \ln(k_B/h) + \Delta S^{\ddagger}/R - (\Delta H^{\ddagger}/R)(1/T)\) and fits by ordinary least squares – the “Eyring plot”: slope \(-\Delta H^{\ddagger}/R\), intercept \(\ln(k_B/h) + \Delta S^{\ddagger}/R\).
- Parameters:
- Return type:
- Returns:
EyringFit
Examples
>>> import numpy as np >>> T = np.linspace(280.0, 360.0, 5) >>> fit = fit_eyring(T, eyring_rate_constant(T, dH=80e3, dS=-20.0)) >>> round(fit.dH, 3), round(fit.dS, 6) (80000.0, -20.0)
- chemistrykit.kinetics.fit_lineweaver_burk(S, v)[source]#
Fit substrate/rate data via the Lineweaver-Burk linearization.
Inverting the Michaelis-Menten equation gives the “double-reciprocal” linear form (Lineweaver & Burk, 1934; Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 20.6):
\[\frac{1}{v} = \frac{K_m}{V_{max}} \cdot \frac{1}{[S]} + \frac{1}{V_{max}}\]so plotting \(1/v\) against \(1/[S]\) gives a straight line of slope \(K_m/V_{max}\) and intercept \(1/V_{max}\), fit here by ordinary least squares.
- Parameters:
- Return type:
- Returns:
MichaelisMentenFit
Examples
Generate exact data from known (Vmax, Km) and recover them:
>>> import numpy as np >>> S = np.array([0.5, 1.0, 2.0, 4.0, 8.0]) >>> v = michaelis_menten_rate(S, Vmax=8.0, Km=1.5) >>> fit = fit_lineweaver_burk(S, v) >>> round(fit.Vmax, 6) 8.0 >>> round(fit.Km, 6) 1.5
- chemistrykit.kinetics.gillespie_ssa(stoich_matrix, rate_constants, reactant_orders, n0, t_max, species=None, seed=None, max_events=10000000)[source]#
Simulate a mass-action network exactly with Gillespie’s direct method.
- Parameters:
stoich_matrix (array-like of int, shape (n_species, n_reactions)) – Net change in each species’ count when each reaction fires.
rate_constants (array-like of float, shape (n_reactions,)) – Stochastic rate constants \(c_j\).
reactant_orders (array-like of int, shape (n_species, n_reactions)) – Number of molecules of each species consumed as reactants by each reaction (0 if not a reactant).
n0 (array-like of int, shape (n_species,)) – Initial molecule counts.
t_max (
float) – Simulate until this time (or until no reaction can fire).species (
Sequence[str] |None) – Species names; defaults to("S0", "S1", ...).seed (int or numpy.random.Generator, optional) – Seed or generator for reproducible trajectories.
max_events (
int) – Safety cap on the number of reaction events.
- Return type:
- Returns:
StochasticTrajectory
Examples
Pure decay
A -> 0of 50 molecules: counts only ever drop by one, and a single trajectory ends with every molecule gone:>>> import numpy as np >>> traj = gillespie_ssa([[-1]], [1.0], [[1]], n0=[50], t_max=100.0, seed=0) >>> int(traj.count("S0")[0]), int(traj.count("S0")[-1]), len(traj.t) (50, 0, 51) >>> bool(np.all(np.diff(traj.count("S0")) == -1)) True
- chemistrykit.kinetics.michaelis_menten_rate(S, Vmax, Km)[source]#
Michaelis-Menten rate law \(v = V_{max}[S] / (K_m + [S])\).
- Parameters:
- Returns:
float or ndarray
Examples
At
S = Kmthe rate is exactly half ofVmax, by construction:>>> round(float(michaelis_menten_rate(S=2.0, Vmax=10.0, Km=2.0)), 6) 5.0
- chemistrykit.kinetics.noncompetitive_inhibition_rate(S, I, Vmax, Km, Ki)[source]#
Noncompetitive-inhibition Michaelis-Menten rate law.
The inhibitor binds a separate site with equal affinity for the free enzyme and the enzyme-substrate complex, which is equivalent to deflating the apparent \(V_{max}\) by \(\alpha = 1 + [I]/K_i\) while leaving \(K_m\) unchanged (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 20.6):
\[v = \frac{V_{max}[S]}{\alpha (K_m + [S])}\]- Parameters:
- Returns:
float or ndarray
Examples
With no inhibitor (
I=0) this reduces exactly to the uninhibited rate law:>>> float(noncompetitive_inhibition_rate(S=2.0, I=0.0, Vmax=10.0, Km=2.0, Ki=1.0)) 5.0
- chemistrykit.kinetics.reversible_analytic(A0, kf, kr, t, B0=0.0)[source]#
Closed-form
A <-> Brelaxation to equilibrium.For \(d[A]/dt = -k_f[A] + k_r[B]\) with \([A]+[B] = A_0+B_0\) conserved, the solution relaxes exponentially to equilibrium with relaxation rate \(k_f + k_r\) (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 20.4):
\[[A](t) = [A]_{eq} + ([A]_0 - [A]_{eq}) e^{-(k_f+k_r)t}, \quad [A]_{eq} = \frac{k_r ([A]_0+[B]_0)}{k_f+k_r}\]- Parameters:
- Returns:
A, B (float or ndarray)
Examples
>>> import numpy as np >>> A, B = reversible_analytic(A0=1.0, kf=2.0, kr=1.0, t=1e6) >>> round(float(A), 3), round(float(B), 3) (0.333, 0.667)
- chemistrykit.kinetics.smoluchowski_rate_constant(D, R_contact)[source]#
Smoluchowski’s steady-state diffusion-limited encounter rate constant.
Two species diffusing with combined diffusion coefficient \(D = D_A + D_B\) react on first contact at separation \(R^* = R_A + R_B\); solving the steady-state diffusion equation around one reactant gives
\[k_D = 4\pi D R^* N_A\]- Parameters:
- Returns:
float – Rate constant, in m^3 mol^-1 s^-1.
Examples
Two small molecules in water (D ~ 2e-9 m^2/s each, contact at 0.5 nm) react at about 1.5e10 L mol^-1 s^-1 (1 m^3 = 1000 L):
>>> k = smoluchowski_rate_constant(D=4e-9, R_contact=5e-10) >>> round(k * 1000 / 1e10, 2) 1.51
- chemistrykit.kinetics.smoluchowski_transient_rate_constant(D, R_contact, t)[source]#
Smoluchowski’s time-dependent diffusion-limited rate constant.
Starting from a uniform (random) distribution of reactants, the depletion zone around each reactant has not yet formed, so the rate coefficient starts higher than \(k_D\) and decays towards it:
\[k(t) = 4\pi D R^* N_A \left(1 + \frac{R^*}{\sqrt{\pi D t}}\right)\]- Parameters:
- Returns:
float or ndarray – Rate constant(s), in m^3 mol^-1 s^-1.
Examples
At \(t = R^{*2}/(\pi D)\) the transient term exactly doubles the steady-state rate:
>>> import numpy as np >>> D, Rc = 4e-9, 5e-10 >>> ratio = smoluchowski_transient_rate_constant(D, Rc, Rc**2 / (np.pi * D)) / smoluchowski_rate_constant(D, Rc) >>> round(float(ratio), 12) 2.0
- chemistrykit.kinetics.ssa_intermediate_concentration(A0, k1, k2, t)[source]#
Steady-state-approximation (SSA) estimate of the intermediate [B].
When the second step is much faster than the first (\(k_2 \gg k_1\)), the intermediate B is consumed almost as fast as it forms, so \(d[B]/dt \approx 0\) at all times after a brief induction period – the steady-state approximation (Espenson, Chemical Kinetics and Reaction Mechanisms, 2nd ed., Ch. 5.1). Applying it to \(d[B]/dt = k_1[A] - k_2[B] = 0\) gives \([B]_{ssa} = (k_1/k_2)[A](t) = (k_1/k_2) A_0 e^{-k_1 t}\), which should agree closely with
consecutive_analytic()’s exact \([B](t)\) whenever \(k_2/k_1\) is large.- Parameters:
- Returns:
float or ndarray
Examples
The SSA error shrinks as k2/k1 grows:
>>> import numpy as np >>> t = np.linspace(0.5, 10.0, 50) >>> errs = [] >>> for k2 in (2.0, 10.0, 100.0): ... _, B_exact, _ = consecutive_analytic(A0=1.0, k1=1.0, k2=k2, t=t) ... B_ssa = ssa_intermediate_concentration(A0=1.0, k1=1.0, k2=k2, t=t) ... errs.append(np.max(np.abs(B_exact - B_ssa))) >>> bool(errs[0] > errs[1] > errs[2]) True