While solving with ode15 i want to put to the limit to variable.

In the code below, Value of (Crec + Ccells) must not exceed Rtot. and after limiting it, the ode must still solve the equations for other t. I need a graph as this.

clc; clear;

opts = odeset('stats','on');

C0 = [6.21E-6, 0, 0, 0, 0]; %Initial Values (umol)
tspan = 0:0.125:4400; %Time span
[t, C] = ode15s(@fn, tspan, C0, opts); %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 on, grid minor
%ylim ([0 7E-6])
%xlim([0 4000])
xlabel('t (in s)'),ylabel('C (in umol/mm^3')
legend('Cmn','Cecm','Crec','Ccirc','Ccells')

function dCdt = fn(t, C)
    %Constant Parameters
    Cmn0 = 6.21E-6; %Initial conc at the MN (umol/MN)
    tr = 1.8E3; %release period (s)
    ka = 5E6; %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 receptors/mm^3)
    r = Cmn0/tr*(t<=tr) + 0*(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

Answer 1

[Alan Stevens]

Like this?

opts = odeset('stats','on');

C0 = [6.21E-6, 0, 0, 0, 0]; %Initial Values (umol)
tspan = 0:0.125:4400; %Time span
[t, C] = ode15s(@fn, tspan, C0, opts); %Solving ODEs in Function (fn

21 successful steps
0 failed attempts
37 function evaluations
1 partial derivatives
6 LU decompositions
27 solutions of linear systems

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 on, grid minor
%ylim ([0 7E-6])
%xlim([0 4000])
xlabel('t (in s)'),ylabel('C (in \mumol/mm^3')
legend('Cmn','Cecm','Crec','Ccirc','Ccells')

function dCdt = fn(t, C)
    %Constant Parameters
    Cmn0 = 6.21E-6; %Initial conc at the MN (umol/MN)
    tr = 1.8E3; %release period (s)
    ka = 5E6; %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 receptors/mm^3)
    r = Cmn0/tr*(t<=tr) + 0*(t>tr);
    
    C_mn  = C(1);
    C_ecm = C(2);
    C_rec = C(3);
    C_circ = C(4);
    C_cells = C(5);
    Csum = C_rec + C_cells;
    if Csum>=Rtot, Csum = Rtot; end
    dCdt = [-r;
             r - ((ka * C_ecm) * (Rtot - Csum)) + (kd * C_rec) - (kc * C_ecm);
             ((ka * C_ecm) * (Rtot - Csum)) - ((kd + ki)* C_rec);
             kc * C_ecm;
             ki * C_rec];
end

[HARSH ZALAVADIYA]

The condition you put, seems correct. Then why isn't I am getting the similar graph to what I have shown? Any idea?

[Alan Stevens]

Possibly you are using different data values.