Simulation of drug distribution with ordinary differential equation

15 次查看（过去 30 天）
Rajvi Amle
Rajvi Amle 2021-5-16
评论： Jeremy Huard 2024-6-12
How can I simulate the code for drug distribution in 10 different random mice at time t=[0 6 12 18 24]. Two variables are used SYSTEM.m_BW and SYSTEM.w_L (range of variables [min max] SYSTEM.m_BW=[24 28], SYSTEM.w_L=[0.419 0.617]). I used the random function in Matlab, which randomly shuffles the mentioned variables. The output function is EXP.Blood_15=(([16.32 3.8 2.87 2.33 2.06]/100).*SYSTEM.m_BW.*SYSTEM.w_L.*DRUG.w_Au_iv.*SYSTEM.m_BW)./(SYSTEM.m_BW. *SYSTEM.w_L). I have also defined the ODE and made a subfile with an ODE equation. I want to introduce a loop in the file toy_ex after the initial condition that calculates EXP.Blood_15 for 10 different iterations/steps (n=10) and save it in a vector y_loop, then to calculate mean and standard deviation from the y_loop. Please help me with this. I have my code here:
function toy_ex
clear all;
close all;
clc;
SYSTEM.m_BW= 24+(28-24).*rand(10,1); % Body weight [g]
SYSTEM.w_L= 0.029+(0.029-0.029).*rand(10,1); % Liver [g/g]
DRUG.w_Au_iv = 15E-03 ; %[mg/g]
EXP.Blood_15=(([16.32 3.8 2.87 2.33 2.06]/100) .*SYSTEM.m_BW.*SYSTEM.w_L.*DRUG.w_Au_iv.*SYSTEM.m_BW)./(SYSTEM.m_BW.*SYSTEM.w_L); %output function
for i=1:n
tspan=[0 6 12 18 24];
MAT=[SYSTEM.m_BW(i) SYSTEM.w_L(i)]
x0=zeros(2,1); % initial condition
[t,x]=ode15s(@ode_ex, tspan, x0, [ ], MAT);
plot(t, x(:,1), 'o', 'Color', blue, 'LineWidth', 2)
end
% ODEs
function dxdt = ode_ex(t, x, SYSTEM)
m_BW=SYSTEM.m_BW;
w_L=SYSTEM.w_L;
m_Au_A=x(1);
m_Au_V=x(2);
% 2 derivatives
dm_Au_A_dt=m_BW - w_L*m_Au_A;
dx(1) = dm_Au_A_dt;
dm_Au_V_dt= w_L - m_BW*m_Au_V;
dx(2) = dm_Au_V_dt;
%vectorized the function
dxdt = dx(:);
end

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2021-5-16
Ran in:
The syntax was not passing the parameters correctly to ‘ode_ex’, and was not passing to correct arguments, at least from what I was able to determine.
This runs without error. See if it does what you want —
SYSTEM.m_BW= 24+(28-24).*rand(10,1); % Body weight [g]
SYSTEM.w_L= 0.029+(0.029-0.029).*rand(10,1); % Liver [g/g]
DRUG.w_Au_iv = 15E-03 ; %[mg/g]
EXP.Blood_15=(([16.32 3.8 2.87 2.33 2.06]/100) .*SYSTEM.m_BW.*SYSTEM.w_L.*DRUG.w_Au_iv.*SYSTEM.m_BW)./(SYSTEM.m_BW.*SYSTEM.w_L); %output function
n = 5;
figure
hold on
for i=1:n
tspan=[0 6 12 18 24];
MAT=[SYSTEM.m_BW(i) SYSTEM.w_L(i)]
x0=zeros(2,1); % initial condition
[t,x]=ode15s(@(t,x)ode_ex(t, x, MAT), tspan, x0);
plot(t, x(:,1), 'o', 'Color', 'blue', 'LineWidth', 2)
end
MAT = 1×2
25.2897 0.0290
MAT = 1×2
27.5717 0.0290
MAT = 1×2
27.6343 0.0290
MAT = 1×2
27.1921 0.0290
MAT = 1×2
26.9200 0.0290
hold off
% ODEs
function dxdt = ode_ex(t, x, SYSTEM)
m_BW=SYSTEM(1);
w_L=SYSTEM(2);
m_Au_A=x(1);
m_Au_V=x(2);
% 2 derivatives
dm_Au_A_dt=m_BW - w_L*m_Au_A;
dx(1) = dm_Au_A_dt;
dm_Au_V_dt= w_L - m_BW*m_Au_V;
dx(2) = dm_Au_V_dt;
%vectorized the function
dxdt = dx(:);
end
.
  6 个评论
显示 4更早的评论
Rajvi Amle
Rajvi Amle 2021-5-21
@Star Strider yes I added ODEs but as I am beginner in Matlab, I started solving easy examples and then add complexities into it step by step.

请先登录，再进行评论。

更多回答（2 个）

Rajvi Amle
Rajvi Amle 2021-5-25
编辑：Rajvi Amle 2021-6-8
Hello @Star Strider. I have a question similarly regarding the code that I asked previously with my first question for simuation of drug distribution with differential equations. In the previous code I had only 2 ODEs but in the present code, I have defined all the variables for the ODEs equations and have now defined 18 ODE equations in a subfile of ODE. Now I want to introduce a loop in the file 'toy_example' after the initial condition that calculates EXP.Blood_15 for 10 different iterations/steps (n=10). I tried it but the code when run shows the error stating that the input arguements are many or less when corrected according to the errors. Also I have created a cell vector to save all the results of x,t in toy_example in a loop to calculate the mean and confidence interval of the simulated data (x(:,2)+x(;,1)). I want to have the graph that runs a loop for 10 iterations (for i=1:10) same as you showed and mentioned on my first question. I worked on it however the code is showing errors for input arguements. Could you please help me with this. I have my code here:
  11 个评论
显示 9更早的评论
Rajvi Amle
Rajvi Amle 2021-6-1
编辑：Rajvi Amle 2021-6-1
Yes. Also, could you please share with me any example for fill function (shaded graphs) for the confidence interval that could be related to my coding regarding the data in a matrix form? Something like given in the following image.
Star Strider
Star Strider 2021-6-1
That is relatively straightforward.
Please provide the data vectors for the regressions and the confidence intervals of the fitted equations.

请先登录，再进行评论。

Rajvi Amle
Rajvi Amle 2021-6-8
Hello @Star Strider. I hope you are doing well. I still have a question regarding confidence interval because I am still confused and not getting significant results. I have 2 dynamic states. x1 and x2 for 10 iterattions or simulations with random variables. Now I have a scenerio with 1 simulation where t1 (time span-n ime steps) is in 6x1 matrix and x (x1 + x2) is in 6X2 matrix. Now I have to calculate the mean for the results x and t from which I can calculate confidence interval. Suppose for matrix x1, the mean should be calculated with the data points given in the row and this average of x1 should be same as time vector t1. And in my original code generated before and also shared with you in the previous questions, I have 25x18 simulation steps which needs to be stored in a matrix. To make you understand better, I have attached a pdf. Please do help me with this. I tried before but not getting through it.
Below is the example of my code to calculate the mean:
function toy_example
clear all;
close all;
clc;
% Body weight and organ weight fraction [Min Max] (random numbers for 10 mice)
%Min=a
%Max=b
%r=a+(b-a)*r(n,1);
n=10;
SYSTEM.m_BW= 24+(28-24).*rand(10,1); % Body weight [g]
SYSTEM.w_L= 0.0417+(0.0681-0.0417).*rand(10,1); % Liver [g/g]
SYSTEM.w_S= 0.019+(0.0051-0.0019).*rand(10,1); % Spleen [g/g]
SYSTEM.w_K= 0.015+(0.0184-0.015).*rand(10,1); % Kidney [g/g]
SYSTEM.w_Lu= 0.0065+(0.0081-0.0065).*rand(10,1); % Lung [g/g]
SYSTEM.w_BR= 0.0139+(0.0191-0.0139).*rand(10,1); % Brain [g/g]
SYSTEM.w_Blood= 0.049+(0.049-0.049).*rand(10,1); % Blood [g/g]
SYSTEM.w_Plasma= 0.029+(0.029-0.029).*rand(10,1); % Plasma [g/g]
%Tissue weight (g)
SYSTEM.m_L = SYSTEM.m_BW.*SYSTEM.w_L; % Liver [g]
SYSTEM.m_BR = SYSTEM.m_BW.*SYSTEM.w_BR; % Brain [g]
SYSTEM.m_K = SYSTEM.m_BW.*SYSTEM.w_K; % Kidney [g]
SYSTEM.m_S = SYSTEM.m_BW.*SYSTEM.w_S; % Spleen [g]
SYSTEM.m_Lu = SYSTEM.m_BW.*SYSTEM.w_Lu; % Lungs [g]
SYSTEM.m_Blood = SYSTEM.m_BW.*SYSTEM.w_Blood; % Blood [g]
SYSTEM.m_Plasma = SYSTEM.m_BW.*SYSTEM.w_Plasma; % Plasma [g]
% Blood flow rate (Fraction of cardiac output)
SYSTEM.F_CC = 16.5; % Cardiac output(mL/h/g) (Upton, 2008; Yuan et al., 2011)
SYSTEM.phi_L = 0.161; % Fraction of blood flow to liver (Average of Buur et al., 2005 and Upton, 2008)
SYSTEM.phi_K = 0.091; % Fraction of blood flow to kidneys (Average of Buur et al., 2005 and Upton, 2008)
SYSTEM.phi_S = 0.011; % Davies and Morris (1993)Table III% Fraction of blood flow to spleen (Average of 6 species)
SYSTEM.phi_BR = 0.033; % Fraction of blood flow to brain (Upton,2008)
%Cardiac output and regional blood flow (mL/h)
SYSTEM.F_C = SYSTEM.F_CC.*SYSTEM.m_BW; % Cardiac output
SYSTEM.F_L = SYSTEM.F_C.*SYSTEM.phi_L; % Blood flow to liver
SYSTEM.F_BR = SYSTEM.F_C.*SYSTEM.phi_BR ; % Blood flow to brain
SYSTEM.F_K = SYSTEM.F_C.*SYSTEM.phi_K ; % Blood flow to kidney
SYSTEM.F_S = SYSTEM.F_C.*SYSTEM.phi_S ; % Blood flow to spleen
SYSTEM.F_Res = SYSTEM.F_C-SYSTEM.F_L-SYSTEM.F_BR-SYSTEM.F_K-SYSTEM.F_S; % Blood flow to rest of body
EXP.t=[0 6 12 18 24]; %h (time for blood pharmacokinetics)
DRUG.w_Au_iv = [15E-03 183E-03 693E-03 1058E-03]; %[mg/g]
DRUG.m_Au_iv = DRUG.w_Au_iv.*SYSTEM.m_BW(1); %[mg]
EXP.Blood_15=(([16.32 3.8 2.87 2.33 2.06]/100).*SYSTEM.m_BW.*SYSTEM.w_Plasma.*DRUG.w_Au_iv(1).*SYSTEM.m_BW)./(SYSTEM.m_BW.*SYSTEM.w_Plasma);
EXP.Blood_183=(([31.3 1.48 0.98 0.84 0.8]/100).*SYSTEM.m_BW.*SYSTEM.w_Plasma.*DRUG.w_Au_iv(2).*SYSTEM.m_BW)./(SYSTEM.m_BW.*SYSTEM.w_Plasma);
EXP.Blood_1058=(([57.82 0.42 0.35 0.35 0.35]/100).*SYSTEM.m_BW.*SYSTEM.w_Plasma.*DRUG.w_Au_iv(4).*SYSTEM.m_BW)./(SYSTEM.m_BW.*SYSTEM.w_Plasma);
t=cell(1,10); y=t; % the output cell arrays as row
figure
hold on
for i=1:10
tspan=[0:24];
MAT=[SYSTEM.m_BW(i) SYSTEM.w_L(i)]';
x0=zeros(2,1);
[t{i},x]=ode15s(@ode_toy, tspan, x0,[], MAT); % assign output to cell arrays
y{i}=(x(:,2)+x(:,1));
plot(t{i},y{i},'-.','Color', 'blue', 'LineWidth', 2)
end
% Confidence interval of 95% for simulated data
mean_y=mean(cell2mat(y),2); % convert cell array to matrix, mean by row
std_y=mean_y/size(y,2); % std dev of mean(y) is mean(y)/nObs
ts = tinv([0.025 0.975],length(y)-1); % T-Score
CI95 = mean_y + ts.*(std_y/sqrt(n));
figure
plot(mean_y) % Plot Mean Of All Experiments
hold on
plot(CI95+mean_y) % Plot 95% Confidence Intervals Of All Experiments
hold off
grid
function dxdt = ode_toy(t,x,SYSTEM)
%% Split State Vector
m_BW=SYSTEM(1);
w_L=SYSTEM(2);
m_Au_A=x(1);
m_Au_V=x(2);
%% ODE system
% Dosing
dm_Au_A_dt=m_BW - w_L*m_Au_A;
dx(1) = dm_Au_A_dt;
dm_Au_V_dt= w_L - m_BW*m_Au_V;
dx(2) = dm_Au_V_dt;
dxdt = dx(:);
end
Please help me with this:
  18 个评论
显示 16更早的评论
Star Strider
Star Strider 2021-6-19
That produces a biphasic plot with an ascending portion and a descending portion. Since you want only the descending portion, it must be isolated from the ascending portion.
This works —
DL1 = load('m_Au_A1.mat');
m_Au_A1 = DL1.m_Au_A1;
DL2 = load('m_Au_V1.mat');
m_Au_V1 = DL2.m_Au_V1;
DL3 = load('t1.mat');
t1 = DL3.t1;
sum1 = m_Au_A1 + m_Au_V1;
[sum1max,idx] = max(sum1);
post_peak = idx:numel(sum1);
t_half_1 = interp1(sum1(post_peak), t1(post_peak), sum1max/2);
figure
plot(t1, sum1)
hold on
plot(t_half_1, sum1max/2, 'sr')
hold off
grid
% set(gca,'XScale','log') % Un-Comment To See The Problem
.
Jeremy Huard
Jeremy Huard 2024-6-12
@Rajvi Amle: I recommend you to have a look into SimBiology which is tailored for (PB)PK modelling. It might solve many of the issues you encounter.
For example it allows to define dosing schedules, provides Monte-Carlo capabilites and percentile plots.
You can use SimBiology both interactively with its Apps or programmatically.
SimBiology Community
Selected users publications with SimBiology and MATLAB (some include SimBiology files)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Biomedical Signal Processing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by