Graphing 1D Diffusion in 2 Compartments - Possible Issue with ODE solver

Hello,

I am using MATLAB to model an experiment where a drug (in a fluid) travels from one compartment (called the "donor compartment"), through a tissue, and into a 2nd compartment (called the "receptor chamber"). You can find a visualization of the system I am modeling below:

I currently have a function called Diffusion_1D_C that models the drug moving out of the donor compartment and through the tissue, and the code can be found below:

function dCdt = Diffusion_1D_C(t, C, x, D1, D2, numX, h1, h, phi)
    dx = x(2) - x(1);
    dCdt = zeros(size(x));
    dCdx = zeros(size(x));
    d2Cdx2 = zeros(size(x));

    % Donor compartment in terms of x coordinates
    nDonor = round((h1 / h) * numX);

    % B.C. zero concentration at bottom of well (i.e., dimensionless x = x/h where x = h, so dimless_x = 1)
    %C(end) = 0;

    % B.C. zero flux (dCdx = 0) at x = 0 - set concentration equal right at top two points, flux = 0
    dCdt(1) = 2 * D1 * (C(2) - C(1))./ (dx.^2);

    % Alpha is right before the interface between donor + tissue
    C_alpha = (D1 * C(nDonor) + (D2 * C(nDonor + 1)))./((D2 * phi) + D1);

    % Beta is right after the interface between donor + tissue
    C_beta = (D1 * C(nDonor) + (D2 * C(nDonor + 1)))./((D1 / phi) + D2);

    for i = 2:nDonor - 1 % 0 to right before h1, needs to be a forwards scheme
        dCdx(i) = (C(i + 1) - C(i))./dx;
        d2Cdx2(i) = (C(i + 1) - 2 * C(i) + C(i - 1))./(dx.^2);
        dCdt(i) = D1 * d2Cdx2(i); % diffusion coeff D1 is for the donor chamber
    end
     
    for i = nDonor % x = last point in donor compartment, use D1, right before interface
        dCdx(i) = (C_alpha - C(i))./ dx;
        d2Cdx2(i) = (C(i - 1) - 2 * C(i) + C_alpha)./(dx.^2);
        dCdt(i) = D1 * d2Cdx2(i);
    end

    for i = nDonor + 1 % x = first point in tissue, use D2, right after interface
        dCdx(i) = (C(i + 1) - C(i))./dx;
        d2Cdx2(i) = ((C_beta) - 2 * C(i) + C(i + 1)) / (dx.^2);
        dCdt(i) = D2 * d2Cdx2(i);
    end

    for i = nDonor + 2:numX - 1 % in tissue compartment, from h1 to h, where h = h1 + h2
        dCdx(i) = (C(i + 1) - C(i)) / dx;
        d2Cdx2(i) = (C(i + 1) - 2 * C(i) + C(i - 1))./(dx.^2);
        dCdt(i) = D2 * d2Cdx2(i);
    end

    for i = numX % very last point, end of tissue, ZERO FLUX BC
        dCdt(i) = D2*((-2*C(i) + C(i - 1))./(dx.^2)); % %BC zero conc
    end 
     
    % make sure dCdt is a column vector..
    dCdt = dCdt(:);
end

In a different function, called varying_thickness_Katz, I graph the concentration in the receptor chamber over time. To do this, I use my Diffusion_1D_C function as well as call ode15s to get the concentration of the drug in the donor compartment and tissue over time. I then find the area under this curve, convert it to mass, and subtract it from the original starting mass of the drug in donor compartment to determine the mass of drug that left the donor compartment and tissue over time. We know that the mass of drug leaving the donor and tissue at a certain time is equal to the mass of drug entering the receptor chamber at a certain time -- we use this to find the mass of drug entering the receptor chamber over time, and subsequently find the concentration of the drug entering the receptor chamber over time. Furthermore, this function also finds confidence intervals around each of the concentration time points in the receptor if there is some technical variability in the thickness of the tissue (labeled as height in the code). I run a certain number of simulations for various tissue heights -- I use normrnd to create these different tissue heights. Once I have these different tissue heights, the concentration in the receptor chamber over time is calculated for a certain number of simulations, and these different concentration profiles are used to find the confidence intervals around each concentration point using a set/defined height. Below is the code:

% Simulation parameters, not dependent on height of tissue 
C0 = 1;  % mol/ml, initial concentration - ratio = 1 at time 0
h1 = 1;  % cm
phi = 1;  % Partition coefficient
tmax = 24 * 60 * 60;  % 24 hours in seconds
numT = 10000;  % Time points
numX = 3001;  % Spatial points
t = linspace(0, tmax, numT);  % Time vector
t_hr = t / 3600;  % Convert to hr
dt = t(2) - t(1);
D1 = 1e-5; % donor compt, low viscosity 
D2 = 1e-5; % tissue
simulations = 10;
tissue_heights = zeros(1,simulations);
All_Conc_receptor = cell(1,simulations);
cIs = zeros(1, simulations); 

%calculate varying thicknesses of tissue, store in array 
for k = 1:simulations
    tissue_heights(k) = normrnd(0.25, 0.05); % in cm, where mean = 0.25 cm, std dev = 0.05 
end 

% create EXPECTED concentration v time plot - h = 0.25 cm 
h_actual = 0.25;
h = h1+h_actual;
nDonor = round((h1 / h) * numX);
x = linspace(0, h, numX);
IC = zeros(size(x));

% initial condition for donor compartment:
IC(1:nDonor) = C0;
% initial condition for tissue compartment
IC(nDonor+1:end) = 0;

% Diffusion model function
diffusion_model = @(t, C) Diffusion_1D_C(t, C, x, D1, D2, numX, h1, h, phi);

% Solve ODE
[t, C] = ode15s(diffusion_model, t, IC);

AUC = zeros(size(1, length(t)));

for i =1:size(C,1)
    % only looking at mass released from donor compt: AUC(i) = trapz(x(1:nDonor),C(i,1:nDonor));
    % look at mass released from donor + tissue compt: 
    AUC(i) = trapz(x(1:end),C(i,1:end));
end 
Mass = AUC.*Area;
M_M0 = Mass./M0;
height_of_donor = 1.0; % cm, change this later.... 
Conc_receptor = (M0 - AUC)/height_of_donor;

% find conc vs time equations for each simulation... 
for i = 1:simulations
    % simulation parameters that depend on height of tissue: 
    h2 = tissue_heights(i);
    h = h2+h1;  % Total distance traveled, cm
    x = linspace(0, h, numX);  % Distance vector
    dx = x(2) - x(1);
    
    M0 = 1;
    Area = pi*(.5)^2;
    C0_M = M0./(Area*h1);

    nDonor = round((h1 / h) * numX);
    IC_M = zeros(size(x));

    % initial condition for donor compartment:
    IC_M(1:nDonor) = C0_M;
    % initial condition for tissue compartment
    IC_M(nDonor+1:end) = 0;
    diffusion_model = @(t, C) Diffusion_1D_C(t, C, x, D1, D2, numX, h1, h, phi);
    [t, C] = ode15s(diffusion_model, t, IC_M);
    % find area under curve 
    AUC = zeros(size(1, length(t)));
    for j =1:size(C,1)
        % find at mass released from donor + tissue compt: 
        AUC(j) = trapz(x(1:end),C(j,1:end));
    end 
    Area_c = h*C0;  % mol/cm^2
    Conc_receptor_temp = (Area_c - AUC)/h;
    All_Conc_receptor{i} = Conc_receptor_temp; 
end 
% find std deviation of concentration at each time point
std_conc = zeros(1, numT);
for sim = 1:simulations
    std_conc = std_conc + (All_Conc_receptor{sim}-mean_concentration).^2;
end
std_conc = sqrt(std_conc/simulations);
%calculation 90% interval
CI = 1.645*std_conc/sqrt(simulations);
% plot mean conc vs. time w/ 90% error bars
figure;
plot(t,Conc_receptor,'b-', 'LineWidth', 2); %plot concentration as func of time
hold on; 
h = errorbar(t, Conc_receptor, CI, 'r', 'LineWidth', 0.25);
% set transparency level
alpha = 0.01;
set([h.Bar, h.Line], 'ColorType', 'truecoloralpha', 'ColorData', [h.Line.ColorData(1:3); 255*alpha])
hold off;

xlabel('Time');
ylabel('Mean Concentration in Receptor Chamber, normalized to 1.0 mol/mL');
title('Mean Concentration vs. Time with 90% Confidence Intervals for Varying Thickness, with Mean Tissue Thickness = 0.25 cm');
legend('Predicted Concentration', 'Location', 'best');                  

When I run the code above using D1 = 1e-5 and D2 = 1e-5, I get results that are positive and make sense according to 1D diffusion in 2 compartments. The graph is shown below:

However, when I run the code using D1 = 1e-6, D2 = 1.e-7 (i.e. make the diffusion coefficients smaller than 1e-5), the concentration plot becomes very weird and does not have the same curve as above, which it should according to 1D diffusion. I think the issue may have to do with the ode solver I am using, since these low of diffusion coefficients may cause it to become inaccurate, but I am not sure which ode solver to use instead OR if the issue comes from somewhere else in my code. Please let me know if you have any advice for how to fix this problem! Also, the graph for D1=1e-6 and D2 = 1e-7 can be found below:

[Torsten]

Is C_W = C_R in your problem formulation ?

You need two transmission conditions to connect C_D and C_R at x=h_D. The condition C_R = phi*C_D is not sufficient to fix a solution.

For "usual" two compartment problems, the transmission conditions would be

C_D = C_R (1) and

D_D*dC_D/dx = D_R*dC_R/dx (2)

(Continuity of concentration and flux at the interface)

So it's not possible only to compute C_D - you will have to compute the spatial and temporal development of C_R as well.

[Charlotte Ellis]

Hey Torsten, thanks so much for the response!

The code I posted has both a concentration and flux boundary conditions at x=h_D. The c_alpha and C_beta concentrations force both the flux to be equal and the partitioning effect of C_R = phi*C_D.

I do compute temporal development of C_R as the code is written.

Any idea on why the code breaks with the second set of diffusion coefficients.

[Torsten]

If you don't compute the temporal development of C_R, it's impossible to force the fluxes to be equal because you don't know the flux of C_R at the interface.

[Charlotte Ellis]

Sorry - I got confused. C_alpha and C_beta ensure flux equality at the boundary between the 2 compartments (donor and tissue) AND the partitioning effect. The downstream boundary condition of the model (between the second cotmpartment and the so called "receptor") is a zero-concentration boundary condition. We compute volume-average concentration in the receptor by first computing mass in the receptor by integrating the flux at the downstream boundary over time then dividing that mass by the volume of the receptor. The receptor is not explicitly in the model.

We compute spacio-temporal concentrations of the donor and tissue compartments by solving the system and compute mass/volume-average concentration in the receptor by integrating the flux at the downstream bounday. Does this make sense?

[Hiep Le]

Can i have 1D DRUG DIFFUSION with the code? The diffusion base on time and space which show effection of the drug on human body

Answer 1

[Torsten]

In my opinion, you have the following conditions at the interface x = h_D:

(1) D_D*dC_D/dx = D_R*dC_R/dx

(2) C_R = phi*C_D (or equivalently dC_R/dt = phi*dC_D/dt assuming C_R(t=0,h_D) = phi*C_D(t=0,h_D) )

(3) dC_D/dt = D_D*d^2C_D/dx^2 (or alternatively dC_R/dt = D_R*d^2C_R/dx^2).

If you introduce ghostpoints xR = h_D + dx_D and xL = h_D - dx_R (xR is right to compartment 1 at distance dx_D and xL is left to compartment 2 at distance dx_R where dx_D,dx_R are the discretization steplengths in compartments 1 and 2, respectively) and discretize equations (1), (2) and (3), you get

(1') D_D*(C_D(xR) - C_D(n-1))/(2*dx_D) = D_R*(C_R^(2)-C_R(xL))/(2*dx_R)

(2') D_R*(C_R(2)-2*C_R(1)+C_R(xL))/(dx_R)^2 = phi * D_D*(C_D(xR)-2*C_D(n)+C_D(n-1))/(dx_D)^2

(3') dC_D(n)/dt = D_D*(C_D(xR)-2*C_D(n)+C_D(n-1))/(dx_D)^2

The equations (1') and (2') determine C_D(xR) and C_R(xL), thus the concentrations in the ghostpoints, and the value for C_D(xR) is used to advance the integration of (3') for C_D(n).

But to compute the gradient of C_R in x = h_d, you must know C_R(2). Thus you must also solve for C_R in the interval h_D <= x <= h_D + h_R.