Biological Effective Dose (BED)

1 view (last 30 days)
Hello, everyone,
Does anyone know how to implement the biological effective dose (BED) calculation? The calculation contains a double integral which I do not know how to calculate.
(#) are number and (v) are vectors.
I am trying to use the cumtrapz.
%% BED
AlphaBeta = 2.6;
nu = 0.3571;
t_w = [0:time_h(end)/length(time_h):time_h(end)]; %??
t_w = t_w';
t_w(end) =[];
Int1 = TAC.*exp(-nu.*(t_w-time_h));
Int2 = TAC.*S.*cumtrapz(time_h,Int1);
G = 2./(max(Dose_Gy3)).^2.*cumtrapz(time_h,Int2);
BED = max(Dose_Gy3).*(1+max(G)./AlphaBeta.*max(Dose_Gy3));
The first internal integral must be made from 0 to t, where t is the time vector. I do not know how to interpret the first integral (dw).
Best regards,
Valentina

Accepted Answer

Torsten
Torsten on 20 Feb 2023
Edited: Torsten on 20 Feb 2023
Assuming you have two arrays (time_h and TAC) of the same length:
AlphaBeta = 2.6;
mu = 0.3571;
S = 1;
int_inner = cumtrapz(time_h,TAC.*exp(mu*time_h));
int_outer = trapz(time_h,TAC.*exp(-mu*time_h).*int_inner);
G = 2./AUC^2*int_outer;
BED = S*AUC.*(1+G./AlphaBeta.*S*AUC);
For example:
format long
mu = 0.3571;;
Ddot = @(t)t.^2;
fun = @(t,omega)Ddot(t).*Ddot(omega).*exp(-mu*(t-omega));
value = integral2(fun,0,1,0,@(t)t)
value =
0.051556597835785
time_h = 0:0.01:1;
TAC = Ddot(time_h);
int_inner = cumtrapz(time_h,TAC.*exp(mu*time_h));
int_outer = trapz(time_h,TAC.*exp(-mu*time_h).*int_inner)
int_outer =
0.051573876925114
  1 Comment
Valentina Vasic
Valentina Vasic on 21 Feb 2023
Dear Torsten,
Perfect, thank you very much.
Exactly what I was looking for. I didn't have in mind splitting the two exponentials. Great idea.

Sign in to comment.

More Answers (1)

Matthieu
Matthieu on 20 Feb 2023
if you have the symbolic formula for TAC :
%Define AUC, S, mu, alpha, beta, T beforehand
syms TAC(t) % Define TAC as symbolic function of t
TAC(t) = t^2 ; % Formula for TAC here
dD_dt(t) = TAC(t)*S ;
D = AUC*S ;
syms t w
int1 = int(dD_dt(w)*exp(-mu*(t-w)),w,0,t) ;
int2 = int(dD_dt(t)*int1,t,0,T) ;
G = (2/D^2)*int2 ;
BED = D*(1+G*D*beta/alpha) ;
use vpa(BED) to get an explicit value for BED.
  1 Comment
Valentina Vasic
Valentina Vasic on 21 Feb 2023
Dear Matthieu,
unfortunately I was not looking for this kind of solution. I have vectors and not defined functions.
But I thank you.

Sign in to comment.

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!