how to use Event function to stop the simulation when certain condition is true with ode15s?

In the code below, the simulation should run with following criteria.

• Crec + Cecm should be smaller than 0.001% Rtot
• Cmn0 and the relative change in Ccells should be less than 0.001%

The above criteria is defined as equilibrim. Therefore the code should run until the equilbrium is achieved. I have checked some answers to similar question earlier, they use 'Event function'. However It is unclear for me, how to use that in my code. Please guide me with it

(NOTE: tspan in not fixed, it can be changed or may be not use if necessary)

clc; clear;

C0 = [6.21E-6, 0, 0, 0, 0]; %Initial Values (umol)
tspan = 0:0.8:4000; %Time span

opts = odeset('Events', @saturation); %Simulation limit
[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
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

And my event function is as following:

function [value,isterminal,direction] = saturation(t,C)

Rtot = 2E-6;

value = (C_rec + C_ecm) - 0.001*Rtot;
insterminal = 1;
direction = 0;

end

Error Message:

Unrecognized function or variable 'C_rec'.

Error in saturation (line 5)
value = (C_rec + C_ecm) - 0.001*Rtot;

Error in odeevents (line 28)
  eventValue = feval(eventFcn,t0,y0,eventArgs{:});

Error in ode15s (line 186)
    odeevents(FcnHandlesUsed,odeFcn,t0,y0,options,varargin);

Error in Original (line 7)
[t, C] = ode15s(@fn, tspan, C0, opts); %Solving ODEs in Function (fn)

Thanks you so much in advance.

Answer 1

[Jan]

Insert these lines:

function [value,isterminal,direction] = saturation(t,C)
Rtot = 2E-6;   % In the other function: Rtot = 1.85E-6 ?!?
C_ecm = C(2);  % <--
C_rec = C(3);  % <--
value = (C_rec + C_ecm) - 0.001*Rtot;
insterminal = 1;
direction = 0;
end

[HARSH ZALAVADIYA]

Hi Jan,

Thank you for replying.

The error message is gone. though it is not the result I was expecting cause I think I have misexplained the question. (Also thank you for point typo in my code)

I have to apply three separate constraints while solving the ode.

1. Crec < 0.001% Rtot
2. Cecm < 0.001% Rtot
3. max(Ccells) < 0.001% Cmn0

Can you guide me how do I attain that?

[Jan]

This is not trivial. You need an event detection to find the time points, where the limits are reached. Then the equations must be changed such, that the derivatives do not exceed 0. As soon as the values are smaller than the limit, the original equations must be enabled again. This is not implemented in Matlab's ODE solvers.