🧬 Math Derivation of CME-defined Stochastic Model of RNA Velocity

Introduction

The stochastic model defined by the Chemical Master Equation (CME) outperforms deterministic ODE models in capturing the inherent stochasticity of single-cell RNA sequencing (scRNA-seq) data. It is actively developed to provide a more accurate representation of feature counts and their underlying biological processes. And it has also enabled the generation of simulated data to evaluate deterministic ODE models and associated data processing methods commonly used in scRNA-seq analysis. Thus I derive the key equations from the paper velocity unraveled, a pivotal paper demonstrating the transformative potential of stochastic approaches.

Overview of Eq(18) from velocity unraveled

🧬 1. Probability generating function (PGF)

We define the state of a cell as:

$$x = (x_u, x_s)$$

Where:

$x_u$ : unspliced mRNA count
$x_s$ : spliced mRNA count

We now write the generating function of their joint distribution:

$$G(u_u, u_s, t) = \sum_x P(x,t)(u_u + 1)^{x_u}(u_s + 1)^{x_s}$$

This is a modified bivariate probability generating function, where the “+1” shift is standard in certain moment-generating setups. It lets you cleanly extract moments via derivatives of $G$ .

⚙️ 2. Characteristic of ODEs derived from CME

$$U_1(u_u, u_s, s) = \frac{u_s \beta}{\beta - \gamma} e^{-\gamma s} + \left(u_u - \frac{u_s \beta}{\beta - \gamma}\right) e^{-\beta s}$$

This expression arises from solving a linear system of ODEs for the chemical master equation (CME) via generating functions, and it essentially encodes how the state propagates in time. It’s derived from how unspliced → spliced reactions occur over time.

🧠 3. Log-generating function $f$

$$f(u_u, u_s, t) := \ln G(u_u, u_s, t) = \int_0^t \alpha(t-s) U_1(u_u, u_s, s) ds$$ This integral form is derived in later section and it tells us how the log of the generating function evolves, driven by transcription rate $\alpha(t-s)$ and the system dynamics encoded in $U_1$ . Essentially, it’s the cumulative effect of production and conversion over time.

Then the log-GF can be written in a linear form in $u_u$ and $u_s$ :

$$f(u_u, u_s, t) = \mu_u(t) u_u + \mu_s(t) u_s$$

This suggests that the process is governed by Poisson distributions, since the log-GF is linear in the arguments.

📊 4. Explicit distribution — product of Poissons

Given the above, we can now recover the joint distribution $P(x,t)$ :

$$P(x,t) = \frac{\mu_u(t)^{x_u} e^{-\mu_u(t)}}{x_u!} \cdot \frac{\mu_s(t)^{x_s} e^{-\mu_s(t)}}{x_s!}$$

So:

$x_u \sim \text{Poisson}(\mu_u(t))$
$x_s \sim \text{Poisson}(\mu_s(t))$
Jointly independent

This model assumes that given time $t$ , the spliced and unspliced counts are independent Poisson-distributed variables, whose rates $\mu_u(t), \mu_s(t)$ evolve in time according to the underlying biochemistry.

🧮 5. Time-averaged distribution

Finally, since single-cell sequencing samples cells asynchronously in time, the observed distribution over counts is not $P(x,t)$ at a fixed $t$ , but a time-averaged version:

$$P(x) = \frac{1}{T} \int_0^T P(x,t) dt$$

This is a mixture of Poissons over time, reflecting the asynchrony of cells in scRNA-seq snapshots. This averaging introduces overdispersion, which is critical to explain the variance observed in real data — greater than what a single Poisson can model.

A) PGF Introduction

1. Multivariate Probability Generating Function (PGF)

This function encodes the entire joint distribution of the random vector $x = (x_u, x_s)$ , i.e., the number of unspliced ( $x_u$ ) and spliced ( $x_s$ ) transcripts.

It’s a bivariate generating function, meaning it’s a function of two complex variables $u_u$ and $u_s$ .

The inclusion of $+1$ makes it a shifted PGF, often done for technical convenience (especially when converting to moment-generating functions).

2. Moment Extraction

From the properties of generating functions:

First Moments:

$$ \frac{\partial G}{\partial u_u}\bigg|_{u_u = u_s = 0} = E[x_u], \quad \frac{\partial G}{\partial u_s}\bigg|_{u_u = u_s = 0} = E[x_s] $$

Second Moments / Covariances:

$$ \frac{\partial^2 G}{\partial u_u^2}\bigg|_{u_u = u_s = 0} = E[x_u(x_u - 1)], \quad \frac{\partial^2 G}{\partial u_u \partial u_s}\bigg|_{u_u = u_s = 0} = E[x_u x_s] $$

These derivatives allow us to compute variances, covariances, and higher-order statistics of transcript counts.

3. Connection to Biological Reactions

In a linear RNA kinetic model:

$x_u$ : produced at rate $\alpha$ , converted to $x_s$ at rate $\beta$
$x_s$ : degraded at rate $\gamma$

The evolution of $G(u_u, u_s, t)$ over time follows a partial differential equation that arises from the Chemical Master Equation (CME), which governs the time evolution of probability distributions in chemical kinetics.

4. Time-Dependence

$G$ is explicitly time-dependent, evolving as the distribution $P(x,t)$ changes.

In some derivations, the log-generating function $f = \log G$ is linear in $u_u, u_s$ , which implies that $x_u, x_s \sim \text{Poisson}(\mu_u(t)), \text{Poisson}(\mu_s(t))$ and are independent.

5. Decay and Stationarity

As $t \to \infty$ :

The means $\mu_u(t), \mu_s(t)$ stabilize.
So does the generating function, converging to that of a product of Poisson distributions (one for each transcript species).

Summary: Why It Matters

Encodes all statistical information about $x_u, x_s$
Enables exact computation of moments, cumulants, and correlations
Links stochastic biochemical kinetics with observed scRNA-seq distributions
Supports modeling of noise, burstiness, and cell-to-cell heterogeneity

B) Derive $u_u(s)$ via method of characteristics

$U_1(u_u, u_s, s)$ is claimed to be $u_u(s)$ , which is the solution to the ODE derived from the Chemical Master Equation (CME) for the two-species birth-death process representing unspliced (u) and spliced (s) mRNA dynamics.

$$U_1(u_u, u_s, s) := u_u(s) = \frac{u_s \beta}{\beta - \gamma} e^{-\gamma s} + \left(u_u - \frac{u_s \beta}{\beta - \gamma}\right) e^{-\beta s}$$

This derivation uses the method of characteristics applied to the generating function of a stochastic process governed by the Chemical Master Equation (CME). We consider a two-species birth-death process representing unspliced (u) and spliced (s) mRNA dynamics.

Step 1: Define the Generating Function

Let the joint probability distribution of unspliced and spliced mRNA at time $t$ be $P(x_u, x_s, t)$ . The generating function is:

$$G(z_u, z_s, t) = \sum_{x_u, x_s} P(x_u, x_s, t) z_u^{x_u} z_s^{x_s}$$

We define new variables:

$$z_u = u_u + 1, \quad z_s = u_s + 1$$

so the generating function becomes:

$$G(u_u, u_s, t) = \sum_{x_u, x_s} P(x_u, x_s, t) (u_u + 1)^{x_u} (u_s + 1)^{x_s}$$

Step 2: CME and Corresponding PDE

The CME for this system is governed by the reactions:

Transcription (birth of unspliced): rate $\alpha$
Splicing: $u \xrightarrow{\beta} s$
Degradation: $s \xrightarrow{\gamma} \emptyset$

From the CME, the PDE for $G$ is below as derived in another section):

$$\frac{\partial G}{\partial t} = \alpha u_u G + \beta(u_s - u_u)\frac{\partial G}{\partial u_u} - \gamma u_s \frac{\partial G}{\partial u_s}$$

Step 3: Method of Characteristics (turn PDE into ODEs)

We now transform this PDE into ODEs using the method of characteristics. Let $u_u(s), u_s(s), G(s)$ be functions of characteristic time $s$ such that along these paths:

$$\frac{du_u}{ds} = \beta(u_s - u_u), \quad \frac{du_s}{ds} = -\gamma u_s, \quad \frac{dG}{ds} = \alpha u_u G$$

Let’s solve these:

Step 4: Solve for $u_s(s)$

$$\frac{du_s}{ds} = -\gamma u_s \Rightarrow u_s(s) = u_s(0) e^{-\gamma s}$$

Step 5: Solve for $u_u(s)$

Use integrating factor method:

$$\frac{du_u}{ds} + \beta u_u = \beta u_s(s) = \beta u_s(0) e^{-\gamma s}$$

Multiply both sides by $e^{\beta s}$ :

$$\frac{d}{ds}(u_u e^{\beta s}) = \beta u_s(0) e^{(\beta - \gamma)s}$$

Integrate:

$$u_u(s) e^{\beta s} = u_u(0) + \frac{\beta u_s(0)}{\beta - \gamma}(e^{(\beta - \gamma)s} - 1)$$

Solve for $u_u(s)$ :

$$u_u(s) = u_u(0) e^{-\beta s} + \frac{\beta u_s(0)}{\beta - \gamma}(e^{-\gamma s} - e^{-\beta s})$$

Step 6: Define $U_1(u_u, u_s, s)$

We identify:

$$U_1(u_u, u_s, s) := u_u(s) = \frac{u_s \beta}{\beta - \gamma} e^{-\gamma s} + \left(u_u - \frac{u_s \beta}{\beta - \gamma}\right) e^{-\beta s}$$

B.1) Derive PDE from CME (step 2 of B)

Step 1: Write the CME explicitly

The CME for the joint distribution $P(x_u, x_s, t)$ of the unspliced $x_u$ and spliced $x_s$ RNA is:

$$ \frac{d}{dt}P(x_u, x_s, t) = \alpha[P(x_u - 1, x_s, t) - P(x_u, x_s, t)] + \beta[(x_u + 1)P(x_u + 1, x_s - 1, t) - x_u P(x_u, x_s, t)] + \gamma[(x_s + 1)P(x_u, x_s + 1, t) - x_s P(x_u, x_s, t)] \qquad (9) $$

where:

$\alpha$ : production rate of unspliced RNA
$\beta$ : splicing rate (unspliced → spliced)
$\gamma$ : degradation rate of spliced RNA

Step 2: Define the generating function $G(u_u, u_s, t)$

$$ G(u_u, u_s, t) := \sum_{x_u=0}^{\infty} \sum_{x_s=0}^{\infty} P(x_u, x_s, t)(u_u + 1)^{x_u}(u_s + 1)^{x_s} \qquad (10) $$

Note: Using $(u_u + 1)^{x_u}$ instead of $u_u^{x_u}$ is a common shift to simplify derivatives later, but you can equivalently define with $u_u^{x_u}$ .

Step 3: Take time derivative of $G$

Using linearity of sums and derivatives:

$$ \frac{\partial G}{\partial t} = \sum_{x_u=0}^{\infty} \sum_{x_s=0}^{\infty} \frac{\partial P(x_u, x_s, t)}{\partial t} (u_u + 1)^{x_u}(u_s + 1)^{x_s} $$

Substitute the CME expression:

$$ = \sum_{x_u,x_s} [\alpha(P(x_u - 1, x_s, t) - P(x_u, x_s, t)) + \beta((x_u + 1)P(x_u + 1, x_s - 1, t) - x_u P(x_u, x_s, t)) + \gamma((x_s + 1)P(x_u, x_s + 1, t) - x_s P(x_u, x_s, t))](u_u + 1)^{x_u}(u_s + 1)^{x_s} $$

Step 4: Evaluate each term separately

Term 1: Transcription

$$ \sum_{x_u,x_s} \alpha(P(x_u - 1, x_s) - P(x_u, x_s))(u_u + 1)^{x_u}(u_s + 1)^{x_s} $$

Rewrite the first sum by shifting $x_u \to x_u + 1$ in the first part:

$$ \sum_{x_u=0}^{\infty} P(x_u - 1, x_s)(u_u + 1)^{x_u} = \sum_{x_u'=-1}^{\infty} P(x_u', x_s)(u_u + 1)^{x_u' + 1} $$

Since $P(x_u', x_s) = 0$ for $x_u' < 0$ , this becomes:

$$ (u_u + 1) \sum_{x_u'=0}^{\infty} P(x_u', x_s)(u_u + 1)^{x_u'} $$

Therefore,

$$ \sum_{x_u,x_s} \alpha(P(x_u - 1, x_s) - P(x_u, x_s))(u_u + 1)^{x_u}(u_s + 1)^{x_s} = \alpha((u_u + 1)G - G) = \alpha u_u G $$

Term 2: Splicing

$$ \sum_{x_u,x_s} \beta((x_u + 1)P(x_u + 1, x_s - 1) - x_u P(x_u, x_s))(u_u + 1)^{x_u}(u_s + 1)^{x_s} $$

Split into two sums:

$$ S_1 = \beta \sum_{x_u,x_s} (x_u + 1)P(x_u + 1, x_s - 1)(u_u + 1)^{x_u}(u_s + 1)^{x_s} $$ $$ S_2 = -\beta \sum_{x_u,x_s} x_u P(x_u, x_s)(u_u + 1)^{x_u}(u_s + 1)^{x_s} $$

Change indices in $S_1$ :

Let $x_u' = x_u + 1$ , $x_s' = x_s - 1$ , so $x_u = x_u' - 1$ , $x_s = x_s' + 1$

Then,

$$ S_1 = \beta \sum_{x_u'=1}^{\infty} \sum_{x_s'=0}^{\infty} x_u' P(x_u', x_s')(u_u + 1)^{x_u' - 1}(u_s + 1)^{x_s' + 1} $$

Rearranged: (see note)

$$ = \beta(u_s + 1) \sum_{x_u',x_s'} x_u' P(x_u', x_s')(u_u + 1)^{x_u' - 1}(u_s + 1)^{x_s'} $$

$S_2$ is:

$$ S_2 = -\beta \sum_{x_u,x_s} x_u P(x_u, x_s)(u_u + 1)^{x_u}(u_s + 1)^{x_s} $$

Now recognize:

$$ \frac{\partial G}{\partial u_u} = \sum_{x_u,x_s} x_u P(x_u, x_s)(u_u + 1)^{x_u - 1}(u_s + 1)^{x_s} $$

So:

$$ S_1 = \beta(u_s + 1) \frac{\partial G}{\partial u_u} $$

and

$$ S_2 = -\beta(u_u + 1) \frac{\partial G}{\partial u_u} $$

Putting together:

$$ \text{Splicing term} = \beta((u_s + 1) - (u_u + 1)) \frac{\partial G}{\partial u_u} = \beta(u_s - u_u) \frac{\partial G}{\partial u_u} $$

Term 3: Degradation

$$ \sum_{x_u,x_s} \gamma((x_s + 1)P(x_u, x_s + 1) - x_s P(x_u, x_s))(u_u + 1)^{x_u}(u_s + 1)^{x_s} $$

Split into:

$$ S_3 = \gamma \sum_{x_u,x_s} (x_s + 1)P(x_u, x_s + 1)(u_u + 1)^{x_u}(u_s + 1)^{x_s} $$ $$ S_4 = -\gamma \sum_{x_u,x_s} x_s P(x_u, x_s)(u_u + 1)^{x_u}(u_s + 1)^{x_s} $$

We have:

$$ S_3 = \gamma \sum_{x_u,x_s} (x_s + 1)P(x_u, x_s + 1)(u_u + 1)^{x_u}(u_s + 1)^{x_s} $$

With the substitution $x_s' = x_s + 1$ , we get $x_s = x_s' - 1$ and:

$$ S_3 = \gamma \sum_{x_u=0}^{\infty} \sum_{x_s'=1}^{\infty} x_s' P(x_u, x_s')(u_u + 1)^{x_u}(u_s + 1)^{x_s' - 1} $$ $$ = \frac{\gamma}{u_s + 1} \sum_{x_u,x_s'} x_s' P(x_u, x_s')(u_u + 1)^{x_u}(u_s + 1)^{x_s'} $$

Recognize

$$ \frac{\partial G}{\partial u_s} = \sum_{x_u,x_s} x_s P(x_u, x_s)(u_u + 1)^{x_u}(u_s + 1)^{x_s - 1} $$ $$ = \frac{\gamma}{u_s + 1} \cdot (u_s + 1) \frac{\partial G}{\partial u_s} = \gamma \frac{\partial G}{\partial u_s} $$

And:

$$ S_4 = -\gamma(u_s + 1) \frac{\partial G}{\partial u_s} $$

So the degradation term is:

$$ \text{Degradation term} = \gamma \frac{\partial G}{\partial u_s} - \gamma(u_s + 1) \frac{\partial G}{\partial u_s} = -\gamma u_s \frac{\partial G}{\partial u_s} $$

Final Result

Combining all terms:

$$ \frac{\partial G}{\partial t} = \alpha u_u G + \beta(u_s - u_u) \frac{\partial G}{\partial u_u} - \gamma u_s \frac{\partial G}{\partial u_s} \qquad (11) $$

B.2) Why Double Sum Equivalence Holds in Term 2 of Step4 in B.1)

This equivalence depends on what values the indices range over and how the function being summed behaves.

1. Notation

When we write:

$$\sum_{x_u', x_s'} f(x_u', x_s')$$

This is shorthand for:

$$\sum_{x_u' = 0}^{\infty} \sum_{x_s' = 0}^{\infty} f(x_u', x_s')$$

That is, summing over all nonnegative integer pairs $(x_u', x_s') \in \mathbb{N}_0 \times \mathbb{N}_0$ .

2. ‘Suspect’ in the derivation

In our original sum:

$$\sum_{x_u' = 1}^{\infty} \sum_{x_s' = 0}^{\infty} x_u' P(x_u', x_s') (u_u + 1)^{x_u' - 1} (u_s + 1)^{x_s'}$$

The lower bound of $x_u'$ is 1 because we performed a change of variables from $x_u = x_u' - 1$ , and in the original sum, $x_u \geq 0$ , which implies $x_u' \geq 1$ .

So the double sum with bounds:

$$\sum_{x_u' = 1}^{\infty} \sum_{x_s' = 0}^{\infty}$$

is not exactly the same as:

$$\sum_{x_u' = 0}^{\infty} \sum_{x_s' = 0}^{\infty}$$

But if we define:

$$\sum_{x_u', x_s'} := \sum_{x_u' = 0}^{\infty} \sum_{x_s' = 0}^{\infty}$$

then in your derivation, the support of $P(x_u', x_s')$ makes this safe because:

$$x_u' P(x_u', x_s') = 0 \text{ when } x_u' = 0$$

since the factor is 0.

So extending the lower limit to 0 adds no contribution to the sum.

✅ Conclusion

The expressions:

$$\sum_{x_u' = 1}^{\infty} \sum_{x_s' = 0}^{\infty} x_u' P(x_u', x_s') \cdots$$

and

$$\sum_{x_u' = 0}^{\infty} \sum_{x_s' = 0}^{\infty} x_u' P(x_u', x_s') \cdots$$

are equal because when $x_u' = 0$ , the term is 0.

Hence, we can write:

$$\sum_{x_u', x_s'} := \sum_{x_u' = 0}^{\infty} \sum_{x_s' = 0}^{\infty}$$

without affecting the value of the sum.

B.3) Solve PDE via Method of Characteristics

We start with the PDE derived from the chemical master equation (CME) for a stochastic model of unspliced (u) and spliced (s) RNA:

$$\frac{\partial G}{\partial t} = \alpha u_u G + \beta(u_s - u_u)\frac{\partial G}{\partial u_u} - \gamma u_s \frac{\partial G}{\partial u_s}$$

This is a first-order linear PDE in 3 variables: $u_u, u_s, t$ .

To solve this, we apply the method of characteristics, which reduces a PDE to a system of ODEs along special curves (characteristics) in the domain $(u_u, u_s, t)$ . The idea is to track how $G$ changes along these curves as we change a parameter $s$ (which can be thought of like an artificial time).

Step 1: Define Characteristic Curves

We introduce $s$ as a parameter along a characteristic curve and define:

$$u_u = u_u(s)$$ $$u_s = u_s(s)$$ $$t = t(s)$$ $$G = G(u_u(s), u_s(s), t(s))$$

Then the total derivative of $G$ along the curve is:

$$\frac{dG}{ds} = \frac{\partial G}{\partial u_u}\frac{du_u}{ds} + \frac{\partial G}{\partial u_s}\frac{du_s}{ds} + \frac{\partial G}{\partial t}\frac{dt}{ds}$$

Now, we substitute the PDE into this expression. From the PDE:

$$\frac{\partial G}{\partial t} = \alpha u_u G + \beta(u_s - u_u)\frac{\partial G}{\partial u_u} - \gamma u_s \frac{\partial G}{\partial u_s}$$

Plugging this into the total derivative:

$$\frac{dG}{ds} = \frac{\partial G}{\partial u_u}\frac{du_u}{ds} + \frac{\partial G}{\partial u_s}\frac{du_s}{ds} + \left[\alpha u_u G + \beta(u_s - u_u)\frac{\partial G}{\partial u_u} - \gamma u_s \frac{\partial G}{\partial u_s}\right]\frac{dt}{ds}$$

Now, choose $\frac{dt}{ds} = 1$ . This simplifies the expression because now $t = s$ , and we can reduce the 3-variable PDE into a system of ODEs in $s$ .

Step 2: Match Terms

To make the right-hand side cancel cleanly, we group terms:

Coefficient of $\frac{\partial G}{\partial u_u}$ :

$$\frac{du_u}{ds} + \beta(u_s - u_u)$$

Coefficient of $\frac{\partial G}{\partial u_s}$ :

$$\frac{du_s}{ds} - \gamma u_s$$

To cancel the dependence on $\frac{\partial G}{\partial u_u}$ and $\frac{\partial G}{\partial u_s}$ , we set these to zero, yielding:

$$\frac{du_u}{ds} = \beta(u_s - u_u), \quad \frac{du_s}{ds} = -\gamma u_s$$

Then the remaining term becomes:

$$\frac{dG}{ds} = \alpha u_u G$$

Now we’ve reduced the PDE into this system of ODEs:

$$\begin{align} \frac{du_s}{ds} &= -\gamma u_s \\ \frac{du_u}{ds} &= \beta(u_s - u_u) \\ \frac{dG}{ds} &= \alpha u_u G \end{align}$$

These are much easier to solve analytically or numerically.

B.4) Legitimacy to choose $\frac{dt}{ds} = 1$ in B.3)

We choose $\frac{dt}{ds} = 1$ in the method of characteristics because it simplifies the partial differential equation (PDE) to a more tractable set of ordinary differential equations (ODEs) — and this choice is completely valid and standard in this method. Let me explain why:

What Does $\frac{dt}{ds} = 1$ Mean

In the method of characteristics, you introduce a new parameter $s$ that traces out a path (or characteristic curve) in the space of independent variables — here, $(u_u, u_s, t)$ . Along this path:

$$\frac{dG}{ds} = \frac{\partial G}{\partial u_u} \frac{du_u}{ds} + \frac{\partial G}{\partial u_s} \frac{du_s}{ds} + \frac{\partial G}{\partial t} \frac{dt}{ds}$$

So $\frac{dt}{ds}$ tells you how “fast” you’re moving in the time direction along the characteristic curve.

Now, if you choose:

$$\frac{dt}{ds} = 1 \Rightarrow t = s$$

you’re saying: “Let the parameter along the path simply equal time.” This simplifies the math without changing the problem.

Why Is This Legitimate

Because in the method of characteristics, $s$ is a dummy variable. You’re free to choose how it relates to the original coordinates, as long as it parametrizes a valid path. The PDE solution is determined by the behavior along these characteristic curves, and the parameterization does not affect the final solution.

Choosing $\frac{dt}{ds} = 1$ :

simplifies the system (fewer variables)
allows you to think of $s$ as time
turns the original PDE into a solvable system of ODEs
does not change the physics or solution of the underlying system

It’s not a constraint imposed by the model — it’s a strategic mathematical choice to simplify the derivation.

Alternative Choices

Yes — you could, in principle, choose something else (e.g. $\frac{dt}{ds} = \gamma$ , or make $s = -t$ ), but then the resulting ODEs are messier. Since the method is agnostic to how you parametrize the path, you’re free to choose the one that leads to the simplest math.

Summary

Choosing $\frac{dt}{ds} = 1$ is valid because:

The method of characteristics allows any parameterization.
This choice makes the equations easier to solve.
The physics or stochastic model (e.g. CME) remains unchanged.

C) Derive log-generating function

Below derive the log-generating function first introduced here.

1. Setup: Use the Method of Characteristics

We rewrite the PDE:

$$\frac{\partial G}{\partial t} = \alpha u_u G + \beta(u_s - u_u)\frac{\partial G}{\partial u_u} - \gamma u_s \frac{\partial G}{\partial u_s}$$

This is a linear PDE in $G(u_u, u_s, t)$ , and we apply the method of characteristics.

Let $s$ be the parameter along characteristic curves. Then we solve the system:

$$\frac{dt}{ds} = 1$$ $$\frac{du_u}{ds} = \beta(u_s - u_u)$$ $$\frac{du_s}{ds} = -\gamma u_s$$ $$\frac{dG}{ds} = \alpha u_u G$$

2. Solve the ODE for $u_s(s)$

$$\frac{du_s}{ds} = -\gamma u_s \Rightarrow u_s(s) = u_s(0)e^{-\gamma s}$$

3. Plug into the ODE for $u_u(s)$

$$\frac{du_u}{ds} = \beta(u_s(s) - u_u(s)) = \beta(u_s(0)e^{-\gamma s} - u_u(s))$$

This is a linear non-homogeneous ODE:

Let’s solve using integrating factor:

Integrating factor: $\mu(s) = e^{\beta s}$

Multiply both sides:

$$e^{\beta s}\frac{du_u}{ds} + \beta e^{\beta s}u_u(s) = \beta u_s(0)e^{(\beta - \gamma)s}$$ $$\frac{d}{ds}(e^{\beta s}u_u(s)) = \beta u_s(0)e^{(\beta - \gamma)s}$$

Integrate both sides:

$$e^{\beta s}u_u(s) = \frac{\beta u_s(0)}{\beta - \gamma}e^{(\beta - \gamma)s} + C$$

Now divide both sides:

$$u_u(s) = \frac{\beta u_s(0)}{\beta - \gamma}e^{-\gamma s} + Ce^{-\beta s}$$

Apply initial condition $u_u(0)$ to solve for $C$ :

$$u_u(0) = \frac{\beta u_s(0)}{\beta - \gamma} + C \Rightarrow C = u_u(0) - \frac{\beta u_s(0)}{\beta - \gamma}$$

So we now have:

$$u_u(s) = \frac{\beta u_s(0)}{\beta - \gamma}e^{-\gamma s} + \left(u_u(0) - \frac{\beta u_s(0)}{\beta - \gamma}\right)e^{-\beta s}$$

4. Define $U_1(u_u, u_s, s)$

Recall that in the characteristic solution for $G$ , we solve:

$$\frac{dG}{ds} = \alpha u_u(s)G \Rightarrow G(s) = \exp\left(\int_0^t \alpha(t-s)u_u(s)ds\right)$$

Define $U_1(u_u, u_s, s) = u_u(s)$ . So:

$$U_1(u_u, u_s, s) = \frac{\beta u_s}{\beta - \gamma}e^{-\gamma s} + \left(u_u - \frac{\beta u_s}{\beta - \gamma}\right)e^{-\beta s}$$

5. Compute $f(u_u, u_s, t) = \ln G$

Now integrate:

$$\ln G(u_u, u_s, t) = \int_0^t \alpha(t-s)U_1(u_u, u_s, s) \, ds$$