Clear Filters
Clear Filters

Biological Effective Dose (BED)

5 views (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.

Categories

Find more on Satellite Mission Analysis in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!