Simulation of drug distribution with ordinary differential equation

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

Answer 1

[Star Strider]

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

.

[Rajvi Amle]

Thank you for your response @Star Strider. Also, how can I calculate and plot the figure for the confidence interval of 95% based on these data if n=10? Could you please help me with this?

[Star Strider]

As always, my pleasure!

Confidence intervals are derived from differences between actual data and fitted regressions, or similarly with respect to known variations in the data. The data here are calculated, with added noise. (It might be preferable to replace rand with randn.) The only way I can think of to add confidence intervals is to somehow combine ‘SYSTEM.m_BW’ and ‘SYSTEM.w_L’ to produce a standard error, with that defined as:

SE = std(x)/sqrt(numel(x));

and use that multiplied by the appropriate result from the tinv function.

I leave those details to you, since I do not understand the calculations.

[Rajvi Amle]

@Star StriderI have a simple question in which I am confused to implement a for loop.I have randomly generated 2 parameters as given above SYSTEM.w_BW and SYSTEM.w_L. Now I have an arbitrary equation as an output function: y=2.34*SYSTEM.w_BW+5.12*SYSTEM.w_L. Now I have to introduce a loop that calculate y for 10 different iterations/steps using (for i=1:10) and save it in a vector y_loop. I also have to calculate mean and std from y_loop. My question is how should I introduce a for loop for different iterations and how to save it in a vector y_loop?

My code is:

clear all;

close all;

clc;

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]

%output function

y=2.34*SYSTEM.m_BW+5.12*SYSTEM.w_L;

for i=1:10;

if i==1

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

if i==2

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

if i==3

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

if i==4

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

if i==5

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

if i==6

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

if i==7

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

if i==8

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

if i==9

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

if i==10

y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);

end

end

I now have to save the loop in a vector y_loop. I am bit confused with this. Could you please help me?

[Star Strider]

First, the if blocks are not necessary, so sliminate them. The addressing in the loop will do this automatically.

The loop then simply becomes —

for i=1:10
    y(i)=2.34*SYSTEM.m_BW(i)+5.12*SYSTEM.w_L(i);
end

Second, if you need to calculate the mean and standard deviation of the ‘y’ vector, simply do that after the end of the loop.

The differential equations and the ode45 call have disappeared. I assume they are no longer necessary?

[Rajvi Amle]

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

[Star Strider]

O.K.

Answer 2

[Rajvi Amle]

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:

[Star Strider]

Unfortunately, that code is so convoluted that I cannot understand it, or what you want to do.

I corrected the ode15s call to:

[t{i},x]=ode15s(@(t,x)ode_toy(t,x,DRUG,SYSTEM), tspan, x0,MAT);        % assign output to cell arrays

however, I cannot figure out the rest of it.

The code now throws a ‘Too many input arguments’ error for that call. The problem is that:

dxdt = ode_toy(t,x,DRUG,SYSTEM)

does not allow for a single ‘MAT’ argument, however that apparenlty is what you want to pass to it.

There is nothing wrong with ‘ode_toy’ that I can see with a quick inspection (not at all thorough, so there could be problems with it), however it may be necessary to create a different function, perhaps something like:

ode_fcn = @(t,x,MAT) ode_toy(t,x,MAT(end-1:end),MAT(1:end-3));

that passes the correct parts of ‘MAT’ to ‘ode_toy’ in the correct order, and then use that with ode15s.

If something like that is what you want to do, then the ode15s call changes to:

[t{i},x]=ode15s(@(t,x)ode_fcn(t,x,MAT), tspan, x0);        % assign output to cell arrays    

This simply illustrates the approach, and is not intended to be definitive code.

I DID NOT TEST ANY OF THIS!

That aside, it would likely be best to begin with a simpler model and then build on it, rather than immediately doing everything at once.

EDIT — Corrected typographical errors.

%% Conversion factors
CF.L_mL=1E3;            %[L] --> [mL]
CF.mg_kg=1E6;           %[mg] --> [kg]
% CF [%ID/g] --> [mg/g]
% 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]
%Organ density (g/cm3,g/mL)- Zhang et al. 2015
SYSTEM.rho_L = 1.06;      % Liver (Average of Buur et al., 2005 and Upton, 2008)
SYSTEM.rho_S = 1.06;      % Spleen (Upton, 2008)
SYSTEM.rho_K = 1.05;      % Kidneys (Average of Buur et al., 2005 and Upton, 2008)
SYSTEM.rho_Lu = 0.26;     % Lungs (Upton, 2008)
SYSTEM.rho_BR = 1.04;     % Brain (Upton, 2008)\
SYSTEM.rho_Blood = 1;     % Blood (Average of Buur et al., 2005 and Upton, 2008)
SYSTEM.rho_Plasma = 1;    % Buur et al. 2005 (Hematocrit is 0.33)
SYSTEM.rho_bile=1;        % Bile
SYSTEM.rho_urine=1;       % Urine
% Tissue volumes (mL) 
SYSTEM.V_L = SYSTEM.m_L./SYSTEM.rho_L;                 % Liver [mL]
SYSTEM.V_BR = SYSTEM.m_BR/SYSTEM.rho_BR;               % Brain [mL]
SYSTEM.V_K = SYSTEM.m_K/SYSTEM.rho_K;                  % Kidney [mL]
SYSTEM.V_S = SYSTEM.m_S/SYSTEM.rho_S;                  % Spleen [mL]
SYSTEM.V_Lu = SYSTEM.m_Lu/SYSTEM.rho_Lu;               % Lungs [mL]
SYSTEM.V_Blood = SYSTEM.m_Blood/SYSTEM.rho_Blood;      % Blood [mL]
SYSTEM.V_Plasma = SYSTEM.m_Plasma/SYSTEM.rho_Plasma;   % Plasma [mL]
% Elimination - Biliary and urinary excretion
SYSTEM.FbileC = 0.00003*CF.L_mL;      %0.0012 [mL/h]
SYSTEM.FurineC = 0.000003*CF.L_mL;    %0.00012[mL/h]
% 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)
% Blood volume fraction in organs and tissues (percentage of tissues)
SYSTEM.phi_BL = 0.31;    % Liver (Buur et al. 2005)
SYSTEM.phi_BS = 0.17;    % Spleen (Brown et al., 1997; Table 30, m_Au_Verage of 3 species)
SYSTEM.phi_BK = 0.24;    % Kidneys (Buur et al., 2005)
SYSTEM.phi_BLu = 0.5;    % Lungs (Brown et al., 1997; Table 30, m_Au_Verage of 3 species)
SYSTEM.phi_BBR = 0.03;   % Brain (Brown et al., 1997; Table 30, m_Au_Verage of 4 species)
SYSTEM.phi_BRes = 0.04;  % Rest of body (Assume the same as the muscle in Buur et al. 2005)
%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
SYSTEM.F_bal = SYSTEM.F_C-SYSTEM.F_L-SYSTEM.F_BR-SYSTEM.F_K-SYSTEM.F_S-SYSTEM.F_Res; %Blood flow balance equation
SYSTEM.F_Bile = SYSTEM.FbileC.*SYSTEM.m_BW;      % [(mL/h)*g]
SYSTEM.F_Urine = SYSTEM.FurineC.*SYSTEM.m_BW;    % [(mL/h)*g]
%
SYSTEM.V_Res = SYSTEM.m_BW-SYSTEM.V_L-SYSTEM.V_BR-SYSTEM.V_K-SYSTEM.V_S-SYSTEM.V_Lu-SYSTEM.V_Plasma;    
SYSTEM.V_bal = SYSTEM.m_BW-SYSTEM.V_L-SYSTEM.V_BR-SYSTEM.V_K-SYSTEM.V_S-SYSTEM.V_Lu-SYSTEM.V_Plasma-SYSTEM.V_Res; 
SYSTEM.V_L_b = SYSTEM.V_L.*SYSTEM.phi_BL;                                                   % Weight/volume of capillary blood in liver compartment
SYSTEM.V_L_t = SYSTEM.V_L-SYSTEM.V_L_b;                                                     % Weight/volume of tissue in liver compartment
SYSTEM.V_BR_b = SYSTEM.V_BR.*SYSTEM.phi_BBR ;                                               % Weight/volume of capillary blood in brain compartment
SYSTEM.V_BR_t = SYSTEM.V_BR-SYSTEM.V_BR_b;                                                  % Weight/volume of tissue in brain compartment
SYSTEM.V_K_b = SYSTEM.V_K.*SYSTEM.phi_BK;                                                   % Weight/volume of capillary blood in kidney compartment
SYSTEM.V_K_t = SYSTEM.V_K-SYSTEM.V_K_b;                                                     % Weight/volume of tissue in kidney compartment
SYSTEM.V_S_b = SYSTEM.V_S.*SYSTEM.phi_BS ;                                                  % Weight/volume of capillary blood in spleen compartment
SYSTEM.V_S_t = SYSTEM.V_S-SYSTEM.V_S_b;                                                     % Weight/volume of tissue in spleen compartment
SYSTEM.V_Lu_b = SYSTEM.V_Lu.*SYSTEM.phi_BLu;                                                % Weight/volume of capillary blood in Lung compartment
SYSTEM.V_Lu_t = SYSTEM.V_Lu-SYSTEM.V_Lu_b;                                                  % Weight/volume of tissue in Lung compartment
SYSTEM.V_Res_b = SYSTEM.V_Res.*SYSTEM.phi_BRes;                                             % Weight/volume of capillary blood in rest of body compartment
SYSTEM.V_Res_t = SYSTEM.V_Res-SYSTEM.V_Res_b;                                               % Weight/volume of tissue in rest of body compartment
%IV dosing
DRUG.t_iv = 0.005;                                       %[h] injection time (approximately 15-20 seconds)
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]
DRUG.M_Au_iv = DRUG.m_Au_iv./DRUG.t_iv;                  %[mg/h]
%Experimental data
EXP.t=[0 6 12 18 24];                                    %h (time for blood pharmacokinetics)
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);
%Enteries in a loop
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) SYSTEM.w_K(i) SYSTEM.w_Lu(i) SYSTEM.w_BR(i) SYSTEM.w_Blood(i) SYSTEM.w_Plasma(i) SYSTEM.m_L(i) SYSTEM.m_BR(i) SYSTEM.m_K(i) SYSTEM.m_S(i) SYSTEM.m_Lu(i) SYSTEM.m_Blood(i) SYSTEM.m_Plasma(i)...
        SYSTEM.V_L(i) SYSTEM.V_BR(i) SYSTEM.V_K(i) SYSTEM.V_S(i) SYSTEM.V_Lu(i) SYSTEM.V_Blood(i) SYSTEM.V_Plasma(i) SYSTEM.F_C(i) SYSTEM.F_L(i) SYSTEM.F_BR(i) SYSTEM.F_K(i) SYSTEM.F_S(i) SYSTEM.F_Res(i) SYSTEM.F_bal(i) SYSTEM.F_Bile(i) SYSTEM.F_Urine(i)...
        SYSTEM.V_Res(i) SYSTEM.V_bal(i) SYSTEM.V_L_b(i) SYSTEM.V_L_t(i) SYSTEM.V_BR_b(i) SYSTEM.V_BR_t(i) SYSTEM.V_K_b(i)  SYSTEM.V_K_t(i) SYSTEM.V_S_b(i) SYSTEM.V_S_t(i) SYSTEM.V_Lu_b(i) SYSTEM.V_Lu_t(i) SYSTEM.V_Res_b(i) SYSTEM.V_Res_t(i)... 
        DRUG.m_Au_iv(1) DRUG.M_Au_iv(1)]';
    options=odeset('AbsTol',10e-2,'RelTol',10e-2,'Stats','on');
    x0=zeros(18,1);
    [t{i},x]=ode15s(@(t,x)ode_toy(t,x,DRUG,SYSTEM), tspan, x0,MAT);        % assign output to cell arrays
    y{i}=(x(:,2)+x(:,1));                           
    plot(t{i},y{i},'-.','Color', 'blue', 'LineWidth', 2)
end

Index exceeds the number of array elements (1).

Error in solution>ode_toy (line 135)
w_L=SYSTEM(2);

Error in solution (line 109)
[t{i},x]=ode15s(@(t,x)ode_toy(t,x,DRUG,SYSTEM), tspan, x0,[]); % assign output to cell arrays

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 152)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

% Confidence interval of 95% for simulated data
n1=1:10;
mean_y=mean(cell2mat(y),2);                 % convert cell array to matrix, mean by row
% N=size(y,1);
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
% nu=N-1;
% conf = 0.95;
% alpha = 1 - conf;
% pLo = alpha/2;
% pUp = 1 - alpha/2;
% crit = tinv([pLo pUp], nu);
CI95 = mean_y + ts.*(std_y/sqrt(n));
figure
plot(mean_y)                                          % Plot Mean Of All Experiments
hold on
plot(n1,CI95+mean_y)                                  % Plot 95% Confidence Intervals Of All Experiments
hold off
grid
function dxdt = ode_toy(t,x,DRUG,SYSTEM)
%% Split State Vector
m_BW=SYSTEM(1); 
w_L=SYSTEM(2); 
w_K=SYSTEM(3);  
w_Lu=SYSTEM(4); 
w_BR=SYSTEM(5);  
w_Blood=SYSTEM(6);  
w_Plasma=SYSTEM(7);  
m_L=SYSTEM(8); 
m_BR=SYSTEM(9);  
m_K=SYSTEM(10);  
m_S=SYSTEM(11); 
m_Lu=SYSTEM(12);  
m_Blood=SYSTEM(13);  
m_Plasma=SYSTEM(14); 
V_L=SYSTEM(15); 
V_BR=SYSTEM(16);  
V_K=SYSTEM(17); 
V_S=SYSTEM(18); 
V_Lu=SYSTEM(19);  
V_Blood=SYSTEM(20); 
V_Plasma=SYSTEM(21);  
F_C=SYSTEM(22);
F_L=SYSTEM(23); 
F_BR=SYSTEM(24);
F_K=SYSTEM(25); 
F_S=SYSTEM(26);
F_Res=SYSTEM(27); 
F_bal=SYSTEM(28); 
F_Bile=SYSTEM(29);  
F_Urine=SYSTEM(30);  
V_Res=SYSTEM(31);  
V_bal=SYSTEM(32); 
V_L_b=SYSTEM(33); 
V_L_t= SYSTEM(34);  
V_BR_b=SYSTEM(35);  
V_BR_t=SYSTEM(36);    
V_K_b=SYSTEM(37);  
V_K_t=SYSTEM(38);  
V_S_b=SYSTEM(39);   
V_S_t=SYSTEM(40);  
V_Lu_b=SYSTEM(41);  
V_Lu_t=SYSTEM(42);  
V_Res_b=SYSTEM(43);  
V_Res_t=SYSTEM(44);
m_Au_iv=DRUG(45);
M_Au_iv=DRUG(46);
m_Au_A=x(1);
m_Au_V=x(2);
m_Au_Lu_b=x(3);
m_Au_Lu_t=x(4);
m_Au_Lu_PC=x(5);
m_Au_BR_b=x(6);
m_Au_BR_t=x(7);
m_Au_Res_b=x(8);
m_Au_Res_t=x(9);
m_Au_K_b=x(10);
m_Au_K_t=x(11);
m_Au_K_PC=x(12);
m_Au_S_b=x(13);
m_Au_S_t=x(14);
m_Au_S_PC=x(15);
m_Au_L_b=x(16);
m_Au_L_t=x(17);
m_Au_L_PC=x(18);
 %% Auxiliary Variables
h=ae_toy(t,DRUG);
%% ODE system
% Dosing
dm_Au_A_dt=F_C*h.rho_Au_V_Lu - F_C*h.rho_Au_A;
dx(1) = dm_Au_A_dt;
dm_Au_V_dt= F_L*(m_Au_L_b./V_L_b) + F_BR*h.rho_Au_V_BR + F_K*h.rho_Au_V_K + F_Res*h.rho_Au_V_Res + h.M_Au_iv - F_C*h.rho_Au_V;
dx(2) = dm_Au_V_dt; 
%% Lung compartment
% Diffusion limited model
dm_Au_Lu_b_dt= F_C*(h.rho_Au_V-h.rho_Au_V_Lu) - h.PALu*h.rho_Au_V_Lu + (h.PALu*h.rho_Au_Lu_t)/PLu;
dx(3) = dm_Au_Lu_b_dt;
dm_Au_Lu_t_dt = h.PALu*h.rho_Au_V_Lu - (h.PALu*h.rho_Au_Lu_t)/PLu-h.r_up_Lu+h.r_rel_Lu;
dx(4) = dm_Au_Lu_t_dt;
dm_Au_Lu_PC_dt = h.r_up_Lu-h.r_rel_Lu;
dx(5) = dm_Au_Lu_PC_dt;
 
%% Brain compartment
% Diffusion limited model
dm_Au_BR_b_dt= F_BR*(h.rho_Au_A-h.rho_Au_V_BR) - h.PABR*h.rho_Au_V_BR + (h.PABR*h.rho_Au_BR_t)/PBR;
dx(6) = dm_Au_BR_b_dt; 
dm_Au_BR_t_dt=h.PABR*h.rho_Au_V_BR - (h.PABR*h.rho_Au_BR_t)/PBR;
dx(7) = dm_Au_BR_t_dt;
%% Rest of body compartment
% Diffusion limited model
dm_Au_Res_b_dt = F_Res*(h.rho_Au_A-h.rho_Au_V_Res) - h.PAres*h.rho_Au_V_Res + (h.PAres*h.rho_Au_Res_t)/Pres;
dx(8) = dm_Au_Res_b_dt; 
dm_Au_Res_t_dt = h.PAres*h.rho_Au_V_Res - (h.PAres*h.rho_Au_Res_t)/Pres;
dx(9) = dm_Au_Res_t_dt;
%% Kidney compartment
% Diffusion limited model
dm_Au_K_b_dt = F_K*(h.rho_Au_A-h.rho_Au_V_K) - h.PAK*h.rho_Au_V_K + (h.PAK*h.rho_Au_K_t)/PK - h.M_Au_Urine;
dx(10) = dm_Au_K_b_dt;
dm_Au_K_t_dt = h.PAK*h.rho_Au_V_K - (h.PAK*h.rho_Au_K_t)/PK-h.r_up_K+h.r_rel_K;
dx(11) = dm_Au_K_t_dt;
dm_Au_K_PC_dt = h.r_up_K-h.r_rel_K;
dx(12) = dm_Au_K_PC_dt;
%% Spleen compartment
% Diffusion limited model
dm_Au_S_b_dt = F_S*(h.rho_Au_A-h.rho_Au_V_S) - h.PAS*h.rho_Au_V_S + (h.PAS*h.rho_Au_S_t)/PS ;
dx(13) = dm_Au_S_b_dt;
dm_Au_S_t_dt = h.PAS*h.rho_Au_V_S - (h.PAS*h.rho_Au_S_t)/PS-h.r_up_S+h.r_rel_S;
dx(14) = dm_Au_S_t_dt;
dm_Au_S_PC_dt = h.r_up_S-h.r_rel_S;
dx(15) = dm_Au_S_PC_dt;
%% Liver compartment
% Diffusion limited model  
% Biliary excretion
dm_Au_L_b_dt = F_L*(h.rho_Au_A-h.rho_Au_V_L) + F_S*h.rho_Au_V_S - h.PAL*h.rho_Au_V_L + (h.PAL*h.rho_Au_L_t)/PL - h.M_Au_Bile ; %- KLRESUP*m_Au_A
dx(16) = dm_Au_L_b_dt;
dm_Au_L_t_dt = h.PAL*h.rho_Au_V_L - (h.PAL*h.rho_Au_L_t)/PL-h.r_up_L+h.r_rel_L;
dx(17) = dm_Au_L_t_dt;
dm_Au_L_PC_dt = h.r_up_L-h.r_rel_L;
dx(18) = dm_Au_L_PC_dt;
dxdt = dx(:);
end
function h=ae_toy(t,x,DRUG,SYSTEM)
%% Split State Vector
m_BW=SYSTEM(1); 
w_L=SYSTEM(2); 
w_K=SYSTEM(3);  
w_Lu=SYSTEM(4); 
w_BR=SYSTEM(5);  
w_Blood=SYSTEM(6);  
w_Plasma=SYSTEM(7);  
m_L=SYSTEM(8); 
m_BR=SYSTEM(9);  
m_K=SYSTEM(10);  
m_S=SYSTEM(11); 
m_Lu=SYSTEM(12);  
m_Blood=SYSTEM(13);  
m_Plasma=SYSTEM(14); 
V_L=SYSTEM(15); 
V_BR=SYSTEM(16);  
V_K=SYSTEM(17); 
V_S=SYSTEM(18); 
V_Lu=SYSTEM(19);  
V_Blood=SYSTEM(20); 
V_Plasma=SYSTEM(21);  
F_C=SYSTEM(22);
F_L=SYSTEM(23); 
F_BR=SYSTEM(24);
F_K=SYSTEM(25); 
F_S=SYSTEM(26);
F_Res=SYSTEM(27); 
F_bal=SYSTEM(28); 
F_Bile=SYSTEM(29);  
F_Urine=SYSTEM(30);  
V_Res=SYSTEM(31);  
V_bal=SYSTEM(32); 
V_L_b=SYSTEM(33); 
V_L_t= SYSTEM(34);  
V_BR_b=SYSTEM(35);  
V_BR_t=SYSTEM(36);    
V_K_b=SYSTEM(37);  
V_K_t=SYSTEM(38);  
V_S_b=SYSTEM(39);   
V_S_t=SYSTEM(40);  
V_Lu_b=SYSTEM(41);  
V_Lu_t=SYSTEM(42);  
V_Res_b=SYSTEM(43);  
V_Res_t=SYSTEM(44);
m_Au_iv=DRUG(45);
M_Au_iv=DRUG(46);
m_Au_A=x(1);
m_Au_V=x(2);
m_Au_Lu_b=x(3);
m_Au_Lu_t=x(4);
m_Au_Lu_PC=x(5);
m_Au_BR_b=x(6);
m_Au_BR_t=x(7);
m_Au_Res_b=x(8);
m_Au_Res_t=x(9);
m_Au_K_b=x(10);
m_Au_K_t=x(11);
m_Au_K_PC=x(12);
m_Au_S_b=x(13);
m_Au_S_t=x(14);
m_Au_S_PC=x(15);
m_Au_L_b=x(16);
m_Au_L_t=x(17);
m_Au_L_PC=x(18);
%% Calculation of Auxiliary Variables
h.M_Au_iv = M_Au_iv*(1.-heaviside(t-DRUG.t_iv));                 %[mg/h]
%Auxiliary variables from ODE
%Blood compartement
h.rho_Au_A = m_Au_A./(V_Plasma*0.2);               %[mg/mL]
h.rho_Au_V=m_Au_V./(V_Plasma*0.8);                 %[mg/mL]
h.rho_Au_V_Lu = m_Au_Lu_b./V_Lu_b;                 %[mg/mL]
h.m_Au_Plasma = m_Au_A+m_Au_V;                      %[mg]
h.w_Au_Plasma=h.m_Au_Plasma./m_Plasma;             %[mg/g]
h.m_Au_Plasmaperc = 100*(h.m_Au_Plasma./m_Au_iv)./(V_Plasma*1000);%[%]
h.m_Au_bloodperc = 100*(h.m_Au_Plasma./m_Au_iv)./(V_Blood*1000);
h.rho_Au_V_L = m_Au_L_b./V_L_b;
h.rho_Au_V_BR = m_Au_BR_b./V_BR_b;
h.rho_Au_V_K = m_Au_K_b./V_K_b;
h.rho_Au_V_Res = m_Au_Res_b./V_Res_b;
h.rho_Au_Plasma = (h.m_Au_Plasma)./(V_Plasma);
%% Lung compartment
% Diffusion limited model
h.rho_Au_Lu_t = m_Au_Lu_t./V_Lu_t;
h.r_up_Lu = h.K_up_Lu*m_Au_Lu_b;
h.r_rel_Lu = K_rel_Lu*m_Au_Lu_PC;
h.m_Au_Lu_total = m_Au_Lu_b+m_Au_Lu_t;
h.rho_Au_Lu = h.m_Au_Lu_total./V_Lu;
h.rho_Au_Lung = (h.m_Au_Lu_total+m_Au_Lu_PC)./V_Lu;
h.rho_Au_Lungtissue = (m_Au_Lu_t+m_Au_Lu_PC)./V_Lu_t;
 
h.m_Au_Lung_Tissue = m_Au_Lu_t+m_Au_Lu_PC;
h.w_Au_Lung_Tissue=h.m_Au_Lung_Tissue./(V_Lu_t*SYSTEM.rho_Lu);
h.w_Au_Lung_total=(m_Au_Lu_t+m_Au_Lu_b+m_Au_Lu_PC)./(V_Lu*SYSTEM.rho_Lu);
h.m_Au_Lungtissueperc = 100*(h.m_Au_Lung_Tissue./m_Au_iv)./(V_Lu_t*1000);
%% Brain compartment
% Diffusion limited model
h.rho_Au_BR_t = m_Au_BR_t./V_BR_t;
h.m_Au_BR_total = m_Au_BR_b+m_Au_BR_t;
h.w_Au_Brain_Tissue=m_Au_BR_t./(V_BR_t*SYSTEM.rho_BR);             %[mg/g]
h.w_Au_Brain_total=(m_Au_BR_t+m_Au_BR_b)./(V_BR*SYSTEM.rho_BR);
h.rho_Au_BR = h.m_Au_BR_total./V_BR;
%% Rest of body compartment
% Diffusion limited model
h.rho_Au_Res_t = m_Au_Res_t./V_Res_t;
h.m_Au_Res_total = m_Au_Res_b+m_Au_Res_t;
h.rho_Au_rest = h.m_Au_Res_total./V_Res;
%% Kidney compartment
% Diffusion limited model
h.r_up_K = h.K_up_K*m_Au_K_b;
h.r_rel_K = K_rel_K*m_Au_K_PC;
h.M_Au_Urine = F_Urine.*h.rho_Au_V_K;      %[mg/h]
h.m_Au_K_total = m_Au_K_b+m_Au_K_t;
h.m_Au_Kidney_Tissue = m_Au_K_t+m_Au_K_PC;
h.rho_Au_K_t = m_Au_K_t./V_K_t;
h.rho_Au_K = h.m_Au_K_total./V_K;
h.rho_Au_Kidney = (h.m_Au_K_total+m_Au_K_PC)./V_K;
h.rho_Au_Kidney_Tissue = (m_Au_K_t+m_Au_K_PC)./V_K_t;
h.w_Au_Kidney_Tissue=h.m_Au_Kidney_Tissue./(V_K_t*SYSTEM.rho_K);
h.w_Au_Kidney_total=(m_Au_K_t+m_Au_K_b+m_Au_K_PC)./(V_K*SYSTEM.rho_K);
h.m_Au_Kidneytissueperc = 100*(h.m_Au_Kidney_Tissue./m_Au_iv)./(V_K_t*1000);
%% Spleen compartment
% Diffusion limited model
h.r_up_S = h.K_up_S*m_Au_S_b;
h.r_rel_S = K_rel_S*m_Au_S_PC;
h.rho_Au_S_t = m_Au_S_t./V_S_t;
h.rho_Au_V_S = m_Au_S_b./V_S_b;
h.m_Au_S_total = m_Au_S_b+m_Au_S_t;
h.rho_Au_S = h.m_Au_S_total./V_S;
h.rho_Au_Spleen = (h.m_Au_S_total+m_Au_S_PC)./V_S;
h.rho_Au_Spleen_Tissue = (m_Au_S_t+m_Au_S_PC)./V_S_t;
h.m_Au_Spleen_Tissue = m_Au_S_t+m_Au_S_PC;
h.w_Au_Spleen_Tissue=h.m_Au_Spleen_Tissue./(V_S_t*SYSTEM.rho_S);
h.w_Au_Spleen_total=(m_Au_S_t+m_Au_S_b+m_Au_S_PC)./(V_S*SYSTEM.rho_S);
h.m_Au_Spleentissueperc = 100*(h.m_Au_Spleen_Tissue./m_Au_iv)./(V_S_t*1000);
%% Liver compartment
% Diffusion limited model 
% Biliary excretion
h.r_up_L = h.K_up_L*m_Au_L_b;
h.r_rel_L = K_rel_L*m_Au_L_PC;
h.M_Au_Bile = F_Bile.*h.rho_Au_V_L; % [mg/h]
h.m_Au_L_total = m_Au_L_b+m_Au_L_t;
h.m_Au_Liver_Tissue = m_Au_L_t+m_Au_L_PC;
h.w_Au_Liver_Tissue=h.m_Au_Liver_Tissue./(V_L_t.*SYSTEM.rho_L);
h.w_Au_Liver_total=(m_Au_L_t+m_Au_L_b+m_Au_L_PC)./(V_L.*SYSTEM.rho_L);
h.rho_Au_L_t = m_Au_L_t./V_L_t;
h.rho_Au_L = h.m_Au_L_total./V_L;
h.rho_Au_Liver = (h.m_Au_L_total+m_Au_L_PC)./V_L;
h.rho_Au_Liver_Tissue = (m_Au_L_t+m_Au_L_PC)./V_L_t;
 
h.m_Au_Livertissueperc = 100*(h.m_Au_Liver_Tissue./m_Au_iv)./(V_L_t.*1000);
% Mass balance
h.Tmass = m_Au_A+m_Au_V+h.m_Au_L_total+h.m_Au_BR_total+h.m_Au_K_total+h.m_Au_Lu_total+h.m_Au_Res_total+h.m_Au_S_total+h.M_Au_Bile+h.M_Au_Urine+m_Au_L_PC+m_Au_S_PC+m_Au_Lu_PC+m_Au_K_PC; %
h.Bal = h.M_Au_iv-h.Tmass;
end

.

[Rajvi Amle]

@Star Strider thank you so much for your response. I understood why the code showed error regarding the input argument. But I did not understand what you said to create a different function 'ode_fcn = @(t,x,MAT) ode_toy(t,x,MAT(end-1:end),MAT(1:end-3));'.

You mean I should create a subfile of 'ode_fcn'?

Also what does this bracket mean "MAT(end-1:end),MAT(1:end-3))"?

[Star Strider]

As always, my pleasure!

My point was to pass a function to ode15s that uses ‘MAT’ as an argument, however that passes elements of ‘MAT’ correctly to to ‘ode_toy’. The current code apparently does not.

[Rajvi Amle]

@Star Strider Hello again.I had a doubt regarding calculation of mean and confidence interval of results x,t. I stored x,t in a cell array to calculate the mean and 95% CI of all the results. But what I did is when I calculated the mean, I defined mean as 'mean_y=mean(cell2mat(y),2);' in the code. Now this states that the mean is calculated for the results of 2nd column instead of both 1st and 2nd columns '(x(:,2)+x(:,1));'. But I want to calculate the mean for the results of both 1st and 2nd columns, like the mean for '(x(:,2)+x(:,1));'.

Following is my code so far:

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
% N=size(y,1);
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

Could you please help me show how to calculate mean of results of both the 1st and 2nd columns of x,t stored in a cell array? That would be really grateful and appreciated.

[Star Strider]

I will do what I can!

This syntax:

mean_y=mean(cell2mat(y),2);

calculates the mean across the columns of ‘y’ for each row. The result is a column vector.

To calculate the column mean values:

mean_y = mean(cell2mat(y));

The column dimension 1 is assumed, so it is not necessary to state it explicitly. This should result in a row vector.

The same applies to the std function that it will be necessary to calculate, in order to calculate the standard error of the mean, and from it the confidence intervals.

std_y = std(cell2mat(y));
std_err_y = std_y/sqrt(size(cell2mat(y),1));

Then use the tinv function to get the critical values of the t-distribution to calculate the confidence intervals for the respective degrees-of-freedom for your data, (since I assume this is how you want to calculate them).

.

[Rajvi Amle]

@Star Strider i tried the code that you showed me

mean_y=mean(cell2mat(y));                           % convert cell array to matrix, mean by row
% N=size(y,1);
std_y = std(cell2mat(y));                           % std dev of mean(y) is mean(y)/nObs
ts = tinv([0.025  0.975],size(y)-1);                % T-Score
CI95 = mean_y + ts.*(std_y./sqrt(size(cell2mat(y),1)));

but still when run, the error stats that

Matrix dimensions must agree.

Error in toy_example
CI95 = mean_y + ts.*(std_y./sqrt(size(cell2mat(y),1)));

I wanted the results to be calculated (mean, std and CI95) in 10x2 matrix stored in a cell array. So how this CI95 can be resolved?

[Star Strider]

I do not understand what is in the (10x2) matrix.

The calculations appear to me to be correct. It may be necessary to loop through them to multiply each of them by ‘ts’,

EDIT — (01 Jun 2021 at 11:40)

Thinking about this further, I believe I now see the problem. The code produced row vectors for the means and standard errors, so that means that the ‘CI95’ calculation should be:

CI95 = mean_y + ts(:)*(std_y./sqrt(size(cell2mat(y),1)));

with ‘ts’ transposed to a column vector with the (:) subscript notation. This produces a (2x10) matrix, so transpose it to get a (10x2) matrix.

[Rajvi Amle]

Thank you @Star Strider. I really appreciate your efforts and time to help me get through these coding problems.

[Star Strider]

As always, my pleasure!

My apologies for not considering that earlier. With row vectors, the ‘EDIT’ approach should work.

[Rajvi Amle]

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]

That is relatively straightforward.

Please provide the data vectors for the regressions and the confidence intervals of the fitted equations.

Answer 3

[Rajvi Amle]

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:

[Rajvi Amle]

Hello @Star Strider. I know it is annoying that I am asking the same question repeatedly. But, as am new on the matlab, I am really confused and have tried a lot doing this code. This is regarding my thesis work so, it would be really very helpful to me if you could suggest me some answers for it. Looking forward to your positive response.

[Star Strider]

I am not certain what other help you want.

Putting this table creation between the code for the figures:

% Output = table(cell2mat(t(1)),mean_y,CI95, 'VariableNames',{'Time','Mean','Confidence_Interval'})

appears to produce the result you want. (I commented it out here so it would not execute in the online Run feature since its arguments have not been defined at this point.)

Other than adding the table, I made a necessary change in the second plot so it would be correct.

toy_example                                 % Call The Function To Run The Code

Output = 25×3 table
    Time     Mean     Confidence_Interval
    ____    ______    ___________________

      0          0          0         0  
      1     25.145     23.346    26.944  
      2     49.021     45.515    52.528  
      3       71.7     66.571     76.83  
      4     93.234     86.564    99.903  
      5     113.68     105.55    121.81  
      6     133.09     123.57    142.61  
      7     151.52     140.68    162.35  
      8     169.01     156.92     181.1  
      9     185.62     172.34     198.9  
     10     201.39     186.99     215.8  
     11     216.37     200.89    231.85  
     12      230.6      214.1    247.09  
     13     244.11     226.65    261.57  
     14     256.94     238.56    275.32  
     15     269.13     249.88    288.39  

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

Output = table(cell2mat(t(1)),mean_y,CI95, 'VariableNames',{'Time','Mean','Confidence_Interval'})


figure
plot(mean_y,'-b')                                     % Plot Mean Of All Experiments
hold on
plot(CI95,'--b')                                      % Plot 95% Confidence Intervals Of All Experiments
hold off
grid
legend('\mu','95% Confidence Limits', 'Location','best')
end


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

.

[Rajvi Amle]

Yes, I wanted the exact graph. But that had only 2 ODEs and the simulation steps were not much. Currently I have results as x = 25x18 simulation steps in my code with 18 ODEs (shown in excel sheet attached). Also,I tried to calculate the confidence interval as suggested by you for the results of x with 10 iterations. But what I found is that the iteration takes place only with the 10th step and does not start from the initial value (the code when run, the workspace shows i=10). As you look into the excel sheet, I want to calculate the mean for two columns as highlighted ((x(:,2)+x(:,1)). But in my code, I think the mean calculated shows only for the first column (also shown in the sheet). Moreover, I am not sure if I did (t=cell(1,10); y=t; % the output cell arrays as row) and my graph of CF for 25x18 simulations as attached in the figure is right. Following is my loop where I am facing problem to plot my iteration steps and the mean.

% t=[];
% x=[];
t=cell(1,10); y=t;                                        % the output cell arrays as row 
figure
hold on

for i=1:n
    tspan=[0:24]; 
    MAT=[SYSTEM.m_BW(i) SYSTEM.w_L(i) SYSTEM.w_K(i) SYSTEM.w_Lu(i) SYSTEM.w_BR(i) SYSTEM.w_Blood(i) SYSTEM.w_Plasma(i) SYSTEM.m_L(i) SYSTEM.m_BR(i) SYSTEM.m_K(i) SYSTEM.m_S(i) SYSTEM.m_Lu(i) SYSTEM.m_Blood(i) SYSTEM.m_Plasma(i)...
        SYSTEM.V_L(i) SYSTEM.V_BR(i) SYSTEM.V_K(i) SYSTEM.V_S(i) SYSTEM.V_Lu(i) SYSTEM.V_Blood(i) SYSTEM.V_Plasma(i) SYSTEM.F_C(i) SYSTEM.F_L(i) SYSTEM.F_BR(i) SYSTEM.F_K(i) SYSTEM.F_S(i) SYSTEM.F_Res(i) SYSTEM.F_bal(i) SYSTEM.F_Bile(i) SYSTEM.F_Urine(i)...
        SYSTEM.V_Res(i) SYSTEM.V_bal(i) SYSTEM.V_L_b(i) SYSTEM.V_L_t(i) SYSTEM.V_BR_b(i) SYSTEM.V_BR_t(i) SYSTEM.V_K_b(i)  SYSTEM.V_K_t(i) SYSTEM.V_S_b(i) SYSTEM.V_S_t(i) SYSTEM.V_Lu_b(i) SYSTEM.V_Lu_t(i) SYSTEM.V_Res_b(i) SYSTEM.V_Res_t(i)... 
        DRUG.m_Au_iv(1) DRUG.M_Au_iv(1)]';
     options=odeset('AbsTol',10e-2,'RelTol',10e-2,'Stats','on');
     x0=zeros(18,1);
     [t{i},x]=ode15s(@ode_toy, tspan, x0,[],ov,DRUG,MAT);     
%         t=[t, {t0}];
%         x=[x, {x0}];
     y{i}=(x(:,2)+x(:,1));                           
     plot(t{i},y{i},'-.','Color', 'blue', 'LineWidth', 2)
end

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

% mean_y=mean(cell2mat(y));                           % convert cell array to matrix, mean by row
% % N=size(y,1);
% std_y = std(cell2mat(y));                           % std dev of mean(y) is mean(y)/nObs
% ts = tinv([0.025  0.975],size(y)-1);                % T-Score
% CI95 = (std_y./sqrt(size(cell2mat(y),1)));

figure
plot(mean_y,'-b')                                          % Plot Mean Of All Experiments
hold on
plot(CI95,'--b')                                    % Plot 95% Confidence Intervals Of All Experiments
hold off
grid
legend('\mu','95% Confidence Limits', 'Location','best')

Could you please just have a look at the loop and the plotting? Suggest me corrections if any.

[Star Strider]

I am lost.

If you want the means of those two columns:

D1 = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/648740/''x''.xlsx', 'Range','D4:E28')

D1 = 25×2

0 0 0.0140 0.0559 0.0100 0.0400 0.0072 0.0290 0.0053 0.0211 0.0039 0.0155 0.0029 0.0114 0.0021 0.0084 0.0016 0.0062 0.0012 0.0046

First2Cols_ColMean = mean(D1(:,[1 2]))

First2Cols_ColMean = 1×2

0.0021 0.0082

First2Cols_RowMean = mean(D1(:,[1 2]),2)

First2Cols_RowMean = 25×1

0 0.0349 0.0250 0.0181 0.0132 0.0097 0.0071 0.0053 0.0039 0.0029

.

[Rajvi Amle]

@Star Strider Is my loop and the plotting in it correct? Because when the code runs, it shows i=10 in the workspace which means only the tenth iteration is stored. Following is my code for loop.

%Enteries in a loop

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) SYSTEM.w_K(i) SYSTEM.w_Lu(i) SYSTEM.w_BR(i) SYSTEM.w_Blood(i) SYSTEM.w_Plasma(i) SYSTEM.m_L(i) SYSTEM.m_BR(i) SYSTEM.m_K(i) SYSTEM.m_S(i) SYSTEM.m_Lu(i) SYSTEM.m_Blood(i) SYSTEM.m_Plasma(i)...
        SYSTEM.V_L(i) SYSTEM.V_BR(i) SYSTEM.V_K(i) SYSTEM.V_S(i) SYSTEM.V_Lu(i) SYSTEM.V_Blood(i) SYSTEM.V_Plasma(i) SYSTEM.F_C(i) SYSTEM.F_L(i) SYSTEM.F_BR(i) SYSTEM.F_K(i) SYSTEM.F_S(i) SYSTEM.F_Res(i) SYSTEM.F_bal(i) SYSTEM.F_Bile(i) SYSTEM.F_Urine(i)...
        SYSTEM.V_Res(i) SYSTEM.V_bal(i) SYSTEM.V_L_b(i) SYSTEM.V_L_t(i) SYSTEM.V_BR_b(i) SYSTEM.V_BR_t(i) SYSTEM.V_K_b(i)  SYSTEM.V_K_t(i) SYSTEM.V_S_b(i) SYSTEM.V_S_t(i) SYSTEM.V_Lu_b(i) SYSTEM.V_Lu_t(i) SYSTEM.V_Res_b(i) SYSTEM.V_Res_t(i)... 
        DRUG.m_Au_iv(1) DRUG.M_Au_iv(1)]';
     options=odeset('AbsTol',10e-2,'RelTol',10e-2,'Stats','on');
     x0=zeros(18,1);
     [t{i},x]=ode15s(@ode_toy, tspan, x0,[],ov,DRUG,MAT);     
     y{i}=(x(:,2)+x(:,1))*(i);                           
     plot(t{i},y{i},'-.','Color', 'blue', 'LineWidth', 2)
end

When I stop the code on the line of plot (as shown above) and debug the code, the iterations starts from i=1,2,3,.... and ends at 10 when the code is continuously made to run. These steps are then stored in y by 25x1 matrix for 10 iteration steps (25x1 steps in a cell of (1,10)) . But, after quitting the debug and the code runs, it shows directly i=10, which means I think only the 10th ietration step is stored in a matrix. Although calculating mean of all the 10 lines (figures in attached files) gives me correct result for confidence interval. Could you please suggest me whether my loop is correct or not? And why only the 10th iteration is stored?

[Star Strider]

I did not run this, however ‘i=10’ indicates that the loop completed successfully.

To determine if all the values were calculated and saved, determine the sizes of ‘t’ and ‘y’. They should have 10 rows and as many columns as necessary.

It may be necessary to use the cell2mat function first to determine the sizes of ‘t’ and ‘y’, and to display them all.

[Rajvi Amle]

The size of t and y when checked in the workplace shows the same size but 1 row and 10 columns - 1x10 cell (and not 10 rows and columns) having 25x1 simulations.

I tried to do t=cell(10,1) and y=t to get 10 rows but then my CF graph is not correct as before.

For cell2mat should I consider another variable like 'a' or 'b' before the loop to determine the sizes of t and y?

[Star Strider]

Using cell2mat will show the the contents of ‘t’ and ‘y’. It may only require veiwing the first and last elements of the cell array to demonstrate that everything exists as it should.

Example —

c{1} = {rand(5)}

c = 1×1 cell array
    {1×1 cell}

m = cell2mat(c{1})

m = 5×5

0.9157 0.3698 0.0598 0.1621 0.8598 0.1127 0.8380 0.5844 0.2586 0.4031 0.5584 0.5496 0.2191 0.1540 0.0166 0.1007 0.4412 0.8340 0.5145 0.2979 0.6890 0.3448 0.3146 0.3537 0.6479

.

[Rajvi Amle]

Thank you for your suggestions @Star Strider. I wanted to ask that how can I calculate the half life based on the above code? is there any differential equation for the half life? As shown in figures of graph attached, I have plotted concentration and area under the curve.Now I want to determine/calculate the half life of given concentration. I have experimental values of half life for the concentration of drug with different doses from the published literature (fig 2: half life.PNG).Can these values be used to plot half life?

[Star Strider]

I do not recall anything that looked like a disappearance curve, so I cannot use something from this example.

In any event, calculating the half-life is straightforward:

t = linspace(0,10);                                 % Create Data
y = 5*exp(-.5*t);                                   % Create Data

t_half = interp1(y, t, max(y)/2)

t_half = 1.3868

figure
plot(t, y)
hold on
plot(t_half, max(y/2), 'rs')
hold off
grid
text(t_half, max(y)/2, sprintf('$\\  \\leftarrow t_{\\frac{1}{2}} = %.2f$', t_half), 'Interpreter','latex')

This is the easiest way, and it does not require any knowledge of the underlying kinetics.

Fitting the data and estimating the parameters and then calculating the half-life from them is also an option.

.

[Rajvi Amle]

@Star Strider I tried plotting half life suggested by you, as data in example from your's was (t,y) , so my data contains (t1,m_Au_V1+m_Au_A1)..... for my code. But in my first graph as attached, the block is not on the curve but inside the curve. I tried managing to bring the block on the curve but didnt work. Could you please have a look on it?

Here is my code so far:

t_half_1 = interp1((m_Au_V1+m_Au_A1),t1, max(m_Au_V1+m_Au_A1)/2);
t_half_2 = interp1((m_Au_V2+m_Au_A2),t2, max(m_Au_V2+m_Au_A2)/2);
t_half_4 = interp1((m_Au_V4+m_Au_A4),t4, max(m_Au_V4+m_Au_A4)/2);

figure(1)
%Blood Pharmacokinetics of 15E-03, 183E-03, 1058E-03 mg/g dose GSH-AuNPs 
set(gcf, 'Position',[0,45,1800,900])
set(0,'defaultAxesFontName', 'Timesnewroman')
set(0,'defaultTextFontName', 'Timesnewroman')
ylabel(['%ID/g in blood(mg/g)'])%CF of (%ID/g) = mg/g
subtightplot(2,2,1, margins)
bar(EXP.t,EXP.Blood_15)
area(t1,(m_Au_V1+m_Au_A1),'FaceColor', [0.9 0.9 0.9])
set(gca,'fontsize',12,'FontName','arial','XScale','lin','Xlim',[0,24],'XTick',[0,6,12,18,24]);
hold on
ylabel(['GSH-AuNPs conc. in blood(mg/g)'])
xlabel(['Time(h)'])
plot(t1,(m_Au_V1+m_Au_A1),'-',EXP.t,EXP.Blood_15,'o',t_half_1,max(m_Au_V1+m_Au_A1)/2,'rs','Color', blue,'LineWidth',2)
er = errorbar(EXP.t,EXP.Blood_15,EXP.errlow_B15,EXP.errhigh_B15);    
er.Color = [0 0 0];                            
er.LineStyle = 'none';
lgd= legend(['AUC = ',num2str(AUC1(end))],'simulated data','experimental data','half-life','Location','northeast');
title(lgd,'15E-03 mg/g ID')
hold off
text(t_half_1, max(m_Au_V1+m_Au_A1)/2, sprintf('$\\  \\leftarrow t_{\\frac{1}{2}} = %.2f$', t_half_1), 'Interpreter','latex')

[Star Strider]

Without having the vectors, I cannot determine what the problem is.

Taking a guess, since the vectors are being added (and I do not know their sizes), try this: —

t_half_1 = interp1((m_Au_V1(:)+m_Au_A1(:)),t1(:), max(m_Au_V1(:)+m_Au_A1(:))/2);

This converts all of them to column vectors, adds them, and then takes half of the maximum of the resulting column, and interpolates to determine the result. (I do not understand why this is not a problem for the others, however they may already be column vectors.)

.

[Rajvi Amle]

Yes sorry @Star Strider. The size of m_Au_A1 is 143x1 and also for t1 and m_Au_V1 it's 143x1. The size for m_Au_A2 is 148x1 and that for m_Au_V2 is 148x1, and its graph for half life is correct but not for (m_Au_A1+m_Au_V1). Also changes with

t_half_1 = interp1((m_Au_V1(:)+m_Au_A1(:)),t1(:), max(m_Au_V1(:)+m_Au_A1(:))/2);

shows no correct graph but same result as before-shows block inside the curve.

[Star Strider]

Please attach (upload) all the vectors either in a .mat file or separate text files. I will need to see them to see what the problem is.

[Star Strider]

I only need the vectors necessary to calculate the half-life. What do I look for?

[Rajvi Amle]

Sorry for the inconvenience @Star Strider. Following are the data for m_Au_A1 and m_Au_V1 and time t1.

[Star Strider]

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]

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