Solving differential equations in Matlab - Modelling Dissolution of a Microsuspension

[poster]

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).

More answer...