MATLAB System of Dependent Differential Equations

Question

Given a system of differential equations such as:

dy/dt = f(t)
dx/dt = g(t)

A solution can be found using dsolve, such as:

dsolve(diff(y) == f(t), diff(x) == g(t), y(0) == 1, x(0) == 1);

But what about a system where all the variables depend on each other:

dy/dt = f(y,z)
dx/dt = g(x,y)
dz/dt = h(z,x)

When approached in the same way, with initial conditions, for a system which does have a solution, I cannot find a solution.

I know the system I have tried can produce solutions as I have used a stochastic/deterministic simulator - think there's probably some strange syntax to use.

I'm specifically looking for the solution where the derivatives are all zero, if that helps.

EDIT:

Here is an example:

PX/dt = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PZ^n))))/kd_mRNA)-kd_prot*PX;
PY/dt = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PX^n))))/kd_mRNA)-kd_prot*PY;
PZ/dt = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PY^n))))/kd_mRNA)-kd_prot*PZ;

with the coefficients:

eff = 20;

KM = 40;

tau_mRNA=2.0;
tau_prot=10;

ps_a=0.5;
ps_0=5.0E-4;

t_ave = tau_mRNA/log(2);

k_tl=eff/t_ave;

a_tr=(ps_a-ps_0)*60;
a0_tr=ps_0*60;

kd_mRNA = log(2)/tau_mRNA;
kd_prot = log(2)/tau_prot;

beta = tau_mRNA/tau_prot; 
alpha = a_tr*eff*tau_prot/(log(2)*KM);
alpha0 = a0_tr*eff*tau_prot/(log(2)*KM);
n=2;

And the initial conditions:

PX0 = 20;
PY0 = 0;
PZ0 = 0;

This produces a response:

This clearly has a steady state solution (all derivatives 0).

In MATLAB I have tried:

%%
syms PX(t) PY(t) PZ(t);

z = dsolve(diff(PX) == (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PZ^n))))/kd_mRNA)-kd_prot*PX, diff(PY) == (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PX^n))))/kd_mRNA)-kd_prot*PY, diff(PZ)==(k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PY^n))))/kd_mRNA)-kd_prot*PZ,PX(0)==20)

and:

%%
eq1 = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PZ^n))))/kd_mRNA)-kd_prot*PX;
eq2 = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PX^n))))/kd_mRNA)-kd_prot*PY;
eq3 = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PY^n))))/kd_mRNA)-kd_prot*PZ;

dsolve(diff(PX)==eq1,PX(0)==20,diff(PY)==eq2,PY(0)==0,diff(PZ)==eq3,PZ(0)==0)

Both produce no errors but return an empty sym.

Answer 1

Your numeric solution appears to have an oscillatory component. The "steady state" may be a zero amplitude limit cycle, which is a non-trivial solution. You definitely shouldn't expect a system like this to have an easy-to-find analytic solution. The cyclic relations between your three variables also doesn't help. For what it's worth, Mathematica 10's DSolve also is unable to find a solution.

Though it won't get you to a solution, the way you're using symbolic math is less than optimal. When you use something like log(2) in a symbolic math equation, 2 should be converted to a symbolic value first. For example, sym(log(2)) yields the approximation 6243314768165359/9007199254740992, whereas log(sym(2)) returns the exact log(2). This latter form is much more likely to lead to solutions if they exist. Here's a modified version of your code, which unfortunately still returns "Warning: Explicit solution could not be found":

eff = 20;

KM = 40;

tau_mRNA=2;
tau_prot=10;

ps_a=1/sym(2);
ps_0=5/sym(10000);

ln2 = log(sym(2));
t_ave = tau_mRNA/ln2;

k_tl=eff/t_ave;

a_tr=(ps_a-ps_0)*60;
a0_tr=ps_0*60;

kd_mRNA = ln2/tau_mRNA;
kd_prot = ln2/tau_prot;

beta = tau_mRNA/tau_prot; 
alpha = a_tr*eff*tau_prot/(ln2*KM);
alpha0 = a0_tr*eff*tau_prot/(ln2*KM);
n=2;

PX0 = 20;
PY0 = 0;
PZ0 = 0;

syms PX(t) PY(t) PZ(t);

eq1 = (k_tl*(a0_tr + a_tr*KM^n/(KM^n + PZ^n))/kd_mRNA)-kd_prot*PX;
eq2 = (k_tl*(a0_tr + a_tr*KM^n/(KM^n + PX^n))/kd_mRNA)-kd_prot*PY;
eq3 = (k_tl*(a0_tr + a_tr*KM^n/(KM^n + PY^n))/kd_mRNA)-kd_prot*PZ;

s = dsolve(diff(PX,t)==eq1,diff(PY,t)==eq2,diff(PZ,t)==eq3,PX(0)==20,PY(0)==0,PZ(0)==0)