The following system of coupled ODEs arises in the study of DNA sequence evolution:
\begin{eqnarray*} \frac{da}{dt} & = & \frac{\mu (1-y) b u}{S - y(S-b-v)} - (\lambda +\mu ) a \\ \frac{db}{dt} & = & -\frac{\mu (b+v) b}{S - y(S-b-v)} + \lambda (1-b) \\ \frac{du}{dt} & = & -\frac{\mu (b+v) u}{S - y(S-b-v)} + \lambda a \\ \frac{dv}{dt} & = & \frac{\mu (b+v) (S-v)}{S - y(S-b-v)} \\ S & = & \exp\left(\frac{\lambda t}{1-x}\right)-1 \end{eqnarray*} for $t>0$, with $a(0)=1$, $b(0)=u(0)=v(0)=0$, $a'(0)=-\lambda-\mu$, $b'(0)=u'(0)=\lambda$, and $v'(0)=0$.
Specifically, they are a moment-matched approximation to the General Geometric Indel model, the simplest continuous-time Markov chain on strings that allows insertions and deletions (a.k.a. "indels") of length > 1; details are here.
The parameters are $\lambda,\mu \in \Re^+$, the indel rates, and $x,y \in [0,1)$, the corresponding probability parameters for the geometric distributions over indel lengths. The mean insertion length is $1/(1-x)$ and the mean deletion length is $1/(1-y)$.
The special case $x=y=0$ (where indels all have unit length) is known as the TKF91 model, and has the following exact solution:
\begin{eqnarray*} a & = & \frac{\mu - \lambda}{\mu L - \lambda M} LM \\ b & = & \frac{\lambda}{\mu L - \lambda M} (L - M) \\ u & = & \frac{\mu - \lambda}{\mu L - \lambda M} M(1-L) \\ v & = & \frac{\mu - \lambda}{\mu L - \lambda M} - 1 \end{eqnarray*} with $L=\exp(-\lambda t)$ and $M=\exp(-\mu t)$.
Can anyone suggest a way to solve this? Things I have tried/observed:
- Clearly $(b,v)$ form an independent subsystem of coupled ODEs, so it's presumably worth solving these first.
- I haven't yet been able to simplify further to just a single ODE.
- I've tried plugging into Mathematica's DSolve as-is. No joy.
- Following the TKF91 solutions, I tried using a trial solution that is a rational function of polynomials that are first-order in $L=\exp(-\lambda t/(1-x))$ and $M=\exp(-\mu t/(1-y))$ and share a common denominator, i.e. $b=\frac{b_0 + b_L L + b_M M + b_{LM} LM}{d_0 + d_L L + d_M M + d_{LM} LM}$ and $v=\frac{v_0 + v_L L + v_M M + v_{LM} LM}{d_0 + d_L L + d_M M + d_{LM} LM}$ (Mathematica gets stuck; if I constrain $d_0=d_{LM}=0$, as in TKF91, it says there is no solution.)
- I can of course approximate it numerically, e.g. with Runge Kutta methods, but I'd really like an analytical solution.
No rush, I've been working on this for 20 years ;)
