Modeling Chemical Reactions using a Stochastic Differential Equation

Question

I am working on a simulation project for a biology professor. Per his advice, I need to use a stochastic differential equation to model the kinetic rate of an organism's metabolism.

Background to the Question:

In a system governed by $N$ biologically catalyzed reactions, here is how I am modeling the reaction rates:

$H_n(t) = X_n(t) \times \prod_{m \in L_n}{C_m}$

Where $H_n$ is the rate of reaction $n$ (in $\mu M$ per day), $C_m$ is the concentration of substrate $m$ (in $\mu M$), and $L_n$ is the set of limiting substrates for reaction $n$. $X_n(t)$ is the kinetic rate coefficient of reaction $n$ at time $t$, which I model with the following Ornstein-Uhlenbeck process:

$X_n(t + h) = \mu + (X_n(t) - \mu)e^{-h\lambda} + f(t)\sigma \sqrt{1 - e^{-2h\lambda}}$

Where $h$ is the step size, $\mu$ is the mean of the stochastic process, $\lambda$ is the rate of decay, $\sigma$ is the standard deviation of the process, and $f(t)$ is Gaussian noise. In my simulation, $t$ and $h$ are in units of days.

In my code, I generate an Ornstein-Uhlenbeck process for all $N$ reactions based on user input. Rather than having the user input values for $\lambda$, $\sigma$, and $\mu$, I have them input values for TRATE (typical reaction rate), TRSTD (typical relative standard deviation), and TDECAY (typical reaction decay).

• TRATE can be loosely interpreted as the mean reaction velocity the rate at which a reaction progresses.
• TDECAY can be interpreted as the rate at which the reaction velocity changes (e.g. the rate of reaction acceleration or deceleration).
• TRSTD is a value between 0 and 1 that represents the % deviation from the typical reaction rate.

For reference, here is a short excerpt from my simulation code:

import numpy as np

# Typical rate of all reactions
TRATE_MIN = 0.01
TRATE_MAX = 1

# Typical decay rate of all Ornstein-Uhlenbeck processes
TDECAY_MIN = 0.01
TDECAY_MAX = 1

# Typical relative standard deviation
TRSTD = 0.1

def calibrate_ornbeck(N) -> list:
    """ Return Ornstein-Uhlenbeck parameters for N processes. """

    assert TRATE_MIN > 0 and TRATE_MAX > 0, "TRATE_MIN and TRATE_MAX must be > 0."
    assert TRATE_MAX > TRATE_MIN, "TRATE_MAX must be strictly greater than TRATE_MIN."

    means = np.zeros(N)
    sigmas = np.zeros(N)
    decays = np.zeros(N)
    starts = np.zeros(N)

    for n in range(N):

        # Values for TRATE pulled from a uniform distribution
        means[n] = np.random.uniform(TRATE_MIN, TRATE_MAX)

        # Values for sigma are based on mean value
        sigmas[n] = np.random.uniform(TRSTD*means[n], means[n]) 

        # Value for TDECAY pulled from a logarithmic distribution
        decays[n] = abs(pow(10, np.random.uniform(np.log10(TDECAY_MIN), np.log10(TDECAY_MAX))))

        # Starting value is picked close to the mean
        starts[n] = abs(np.random.normal(means[n], sigmas[n]))

    # return mu, sigma, lambda, and starting values (X_n(0))
    return means, sigmas, decays, starts

After having the user input values for TRATE, TDECAY, and TRSTD I then generate values for the Ornstein-Uhlenbeck processes based on the parameters returned from calibrate_ornbeck. When I execute the simulation, I interpolate between those values using a spline function.

My Question:

Currently, I have set TRATE and TDECAY to be some value between 0.01 and 1. However, I want to expand this range of values to account for variability in the kinetic rates. I am not sure what I should set them to though, and I'm having a hard time finding some good papers on how I can set these values. I've found a lot of papers on how OU processes have been used to model changes in gene expression, but nothing closely resembling the kind of analysis my professor would like me to do.

Would anyone have any advice about how to model chemical reactions using a stochastic process? I feel my mentor has helped me determine the right mathematical framework, but I am lost as to how to reasonably parameterize it.