Graphing 1D Diffusion in 2 Compartments - Possible Issue with ODE solver
2 次查看(过去 30 天)
显示 更早的评论
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1660196/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1660201/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1660206/image.png)
5 个评论
回答(1 个)
Torsten
2024-4-11
编辑:Torsten
2024-4-11
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.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!