Solving differential equations in Matlab

Question

I need to solve these 2 differential equations simultaneously.

dr^3/dt=(-3*D*Cs)/(ρ*r0^2 )*r*(1-C)

dC/dt=((D*4π*r0*N*(1-C)*r)-(Af*C))/V

Note: dr^3/dt is the derivative of r^3 with respect to t

The two equations resemble the change in particle radius (r) and concentration (C) with time for a dissolution process of a microsuspension and its simultaneous absorption in the bloodstream. What is expected to happen as the solid dissolves, is that radius, r, will decrease and the concentration, C, will increase and eventually plateau (i.e. reach an equilibrium) as the dissolved solid is being removed into the bloodstream by this Af*C term (where Af is some sort of absorption rate constant). The equations come from this paper which I am trying to replicate: jpharmsci.org/article/S0022-3549(18)30334-4/fulltext#sec3.2.1 -- Change in C with t is supposed to be like Figure 3 (DCU example).

I did the simplification: dr^3/dt = 3r^2*(dr/dt) and by dividing both sides of the equation by 3r^2. The odes become:

function dydt=odefcnNY_v3(t,y,D,Cs,rho,r0,N,V,Af)
dydt=zeros(2,1);
dydt(1)=((-D*Cs)/(rho*r0^2*y(1)))*(1-y(2)); % dr*/dt
dydt(2)=(D*4*pi*N*r0*(1-y(2))*y(1)-(Af*y(2)))/V; % dC*/dt
end

y(1) = r* and 
y(2) = C*

r* and C* 

is the terminology used in the paper and are "normalised" radius and concentration where

r*=r/r0 and C*=C/Cs

where:

• r=particle radius (time varying and expressed by dydt(1))
• r0=initial particle radius
• C=concentration of dissolved solids (time varying and expressed by dydt(2))
• Cs=saturated solubility

The rest of the code is below. Updated with feedback from authors on values used in paper and to correct initial values to y0=[1 0]

MW=224; % molecular weight
D=9.916e-5*(MW^-0.4569)*60/600000 %m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60] equation provided by authors, divide by 600,000 to convert to m2/s 
rho=1300; %kg/m3
r0=10.1e-6; %m dv50
Cs=1.6*1e6/1e9; %kg/m3 - 1.6ug/m3 converted to kg/m3
V=5*0.3/1e6;%m3 5ml/Kg animal * 0.3Kg animal, divide by 1e6 to convert to m3
W=30*0.3/1000000; %kg; 30mg/Kg animal * 0.3Kg animal, divide by 1e6 to convert to m3
N=W/((4/3)*pi*r0^3*rho); % particle number
Af=0.7e-6/60; %m3/s
tspan=[0 24*3600]; %s in 24 hrs
y0=[1 0];
[t,y]=ode113(@(t,y) odefcnNY_v11(t,y,D,Cs,rho,r0,Af,N,V), tspan, y0);
plot(t/3600,y(:,1),'-o') %plot time in hr, and r*
xlabel('time, hr')
ylabel('r*, (rp/r0)')
legend('DCU')
title ('r*');
plot(t/3600,y(:,1)*r0*1e6); %plot r in microns
xlabel('time, hr');
ylabel('r, microns');
legend('DCU');
title('r');
plot(t/3600,y(:,2),'-') %plot time in hr, and C*
xlabel('time, hr')
ylabel('C* (C/Cs)')
legend('DCU')
title('C*');
plot(t/3600, y(:,2)*Cs) % time in hr, and bulk concentration on y
xlabel('time, hr')
ylabel('C, kg/m3')
legend('Dissolved drug concentration')
title ('C');

I first tried ode45, but the code was taking a very long time to run and eventually I got some errors. I then tried ode113 and got the below error.

Warning: Failure at t=2.112013e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

Update: Code for function updated to resolve singularity issue:

function dydt=odefcnNY_v10(t,y,D,Cs,rho,r0,N,V,Af)
dydt=zeros(2,1);
dydt(1)=(-D*Cs)/(rho*r0^2)*(1-y(2))*y(1)/(1e-6+y(1)^2); % dr*/dt
dydt(2)=(D*4*pi*N*r0*(1-y(2))*y(1)-Af*y(2))/V; %dC*/dt
end

Results

Mechanistic background to model

Derivation of dr/dt

Derivation of dC/dt

Model Assumptions

Above you will find slides that show the derivations of these equations. They assumed that in the Noyes-Whitney equation for dissolution rate dM/dt=(/h)4^2(−), the film thickness, h, is equal to the particle radius, r. This is a simplification usually done in biopharmaceutics modelling if the Reynolds number is low and particles are <60um (in their case 10um). If we make this assumption, we are left with dM/dt=4(−). I am eager to replicate this paper as I want to do this exact same thing i.e. model drug absorption of a subcutaneous injection of a microsuspension. I have contacted the authors, who seem rather unsure of what they have done so I am looking at other sources, there is for example this paper: https://pubs.acs.org/doi/pdf/10.1021/acs.iecr.7b04730 where in equation 6, an equation for dC/dt is shown. They imbed the change in surface area per unit volume (a) (equation 5) into equation 6. And their mass transfer coefficient kL is a lumped parameter = D/h (diffusivity/film thickness).

Answer 1

In the original form of the radius equation

d(r^3)/dt = -3K*(r^3)^(1/3)*(1-C)

or the power-reduced one

dr/dt = -K/r*(1-C)  <==> d(r^2)/dt = -2K*(1-C)

you reach a singularity at the moment that the radius shrinks to zero like approximately r=sqrt(A-B*t), that is, the globules will have vanished in finite time. From that point on, obviously, one should have r=0. This can be obtained via a model change from a 2 component system to a scalar equation for C alone, getting the exact time via event mechanism.

Or alternatively, you can modify the radius equation so that r=0 is a natural stationary point. The first version somehow does that, but as the second versions are equivalent, one can still expect numerical difficulties. Possible modifications are

d(r^2)/dt = -2K*sign(r)*(1-C)

where the signum function can be replaced by continuous approximations like

x/(eps+abs(x)), x/max(eps,abs(x)), tanh(x/eps), ...

or in the reduced form, the singularity can be mollified as

dr/dt = -K*(1-C) * r/(eps^2+r^2)

or still in some other variation

dr/dt = -K*(1-C) * 2*r/(max(eps^2,r^2)+r^2)

---

Applied to the concrete case this gives (using python instead of matlab)

def odefcnNY(t,y,D,Cs,rho,r0,N,V,Af):
    r,C = y;
    drdt = (-D*Cs)/(rho*r0**2)*(1-C) * r/(1e-6+r**2); # dr*/dt
    dCdt = (D*4*pi*N*r0*(1-C)*r-(Af*C))/V;            # dC*/dt
    return [ drdt, dCdt ];

and apply an implicit method in view of the near singularity at r=0

D=8.3658e-10#m2/s 
rho=1300; #kg/m3
r0=10.1e-6; #m dv50
Cs=0.0016; #kg/m3
V=1.5e-6;#m3
W=9e-6; #kg
N=W/(4/3*pi*r0^3*rho);
Af=0.7e-6/60; #m3/s
tspan=[0, 24*3600]; #sec in 24 hours
y0=[1.0, 0.0]; # relative radius starts at full, 1.0
sol=solve_ivp(lambda t,y: odefcnNY(t,y,D,Cs,rho,r0,Af,N,V), tspan, y0, method="Radau", atol=1e-14);
t = sol.t; r,C = sol.y;

then results in a solution plot

which now with the corrected parameters looks close to the published graphs.