Solving differential equations to get accurate plot as shown.

I have five odes with initial conditions to plot. I am anticipating the graph as shown.

There is also condition of (Crec + Ccells) must not exceed (Rtot). I dont know how to put this. My code is below.

clc; clear; close;

C0 = [6.2E-6, 0, 0, 0, 0]; %Initial Values
tspan = 0:1:10000; %Time span
[t, C] = ode23t(@fn, tspan, C0); %Solving ODEs in Function (fn)

C_mn = C(:,1);
C_ecm = C(:,2);
C_rec = C(:,3);
C_circ = C(:,4);
C_cells = C(:,5);

plot(t,C_mn,t,C_ecm,t,C_rec,t,C_circ,t,C_cells),grid
xlim([0 10000])
xlabel('t'),ylabel('C')
legend('Cmn','Cecm','Crec','Ccirc','Ccells')

function dCdt = fn(t, C)
    %Constant Parameters
    Cmn0 = 6.2E-6; %Initial conc at the MN (
    tr = (3600+10)/2; %release period (s)
    ka = 5.01E6; %Association rate (in 1/s)
    kd = 5E-4; %Dissociation rate (in 1/s)
    ki = 5.05E-3; %Internalization rate (in 1/s)
    kc = 5E-3; %Circulation uptake rate (in 1/s)
    Rtot = 1.85E-6 ; %initial receptor concentration (in umol/mm^3)
    r = Cmn0/tr*(t<=tr);
    
    C_mn  = C(1);
    C_ecm = C(2);
    C_rec = C(3);
    C_circ = C(4);
    C_cells = C(5);
    
    dCdt = [-r;
             r - ((ka * C_ecm) * (Rtot - C_rec - C_cells)) + (kd * C_rec) - (kc * C_ecm);
             ((ka * C_ecm) * (Rtot - C_rec - C_cells)) - ((kd + ki)* C_rec);
             kc * C_ecm;
             ki * C_rec];
end

I am getting this

Instead of

[VBBV]

Use

ylim([0 8e-6])

[VBBV]

tspan = 0:1:3000;

Use tspan to solve ode

[VBBV]

C0 = [1.6E-6, 0, 0, 0, 0];

Use a good initial guess

Answer 1

[Bjorn Gustavsson]

Since you get an oscillatory behaviour of (primarily) Cecm and Crec it is very likely that you have a bug for some of the reactions. You might also strongly benefit from using the NonNegative option forcing the solution to be positive (I assume negative concentrations are as unchemical as they are unphysical), to do that have a look at the help and documentation of odeset.

HTH