# funzioni per calcoli "biologici"...
# VARIABILI DI "SISTEMA" -------------------------------------------------------
#' srim.limits
srim.let.max <- c(H=8.431E-02, He=2.369E-01, Li=4.040E-01, Be=6.052E-01, B=7.825E-01, C=9.757E-01, N=1.159E+00, O=1.376E+00, F=1.609e+00, Ne=1.689e+00)*1000
#' srim.limits
srim.let.min <- c(H=2.060E-04, He=8.243E-04, Li=1.900E-03, Be=3.351E-03, B=5.484E-03, C=8.126E-03, N=1.174E-02, O=1.606E-02, F=1.964e-02, Ne=2.458e-02)*1000
#' srim.limits
srim.es.max <- c(H=10.00E+03, He=2.500e+03, Li=1.428e+03, Be=1.111e+03, B=9.090e+02, C=8.333e+02, N=7.142e+02, O=6.2500e+02, F=6.842680e+02, Ne=6.442118e+02)
#' srim.limits
srim.es.min <- c(H=10.00E-03, He=2.500e-03, Li=1.428e-03, Be=1.111e-03, B=9.090e-04, C=8.333e-04, N=7.142e-04, O=6.2500e-04, F=5.263600e-04, Ne=4.955475e-04)
#' MKM lower limits
mkm.lower <- c(0.1,0.02,1,0.1)
#' MKM upper limits
mkm.upper <- c(1,1,10,1)
# BASIC BIOLOGICAL MODELS ------------------------------------------------------
#' Evaluation of Linear-Quadratic (LQ) Poisson Model parameters
#'
#' It translates R, D50 and gamma parameters to alphaX and betaX.
#'
#' Note: d.ref è la dose per frazione di riferimento
#' in cui sono fittati i dati per D50, gam.
#' Per il momento d.ref non è usato: si assume una condizione di riferimento
#' standard per D50 e gam (forse... )
#'
#' @param R the ratio alphaX/betaX (Gy).
#' @param D50 the tolerance dose at 50\%.
#' @param gam the gamma parameter (associated to the clonogenic density).
#' @param d.ref the reference dose at which D50 and gam are fitted (usually 2 Gy).
#' @param d the actual level of dose for the evaluation of alphaX, betaX.
#'
#' @return A list consisting of:
#' \item{alphaX}{the LQ alpha parameter}
#' \item{betaX}{the LQ beta parameter}
#'
#' @family Basic Biological Models
#' @export
alpha.beta.lqp <- function(R=R, D50=D50, gam=gam, d.ref=2, d=2)
{
e <- exp(1)
D50 <- D50*(1+d.ref/R)/(1+d/R)
alphaX <- (e*gam-log(log(2)))/(D50*(1+d/R))
betaX <- (e*gam-log(log(2)))/(D50*(R+d))
return(list(alphaX=alphaX, betaX=betaX))
}
#' Evaluation of dose-depending RBE
#'
#' Evaluate the Relative Biological effectiveness (RBE) of the actual radiation.
#' It uses the parameters: alpha, beta, alphaX and betaX
#'
#' @param alpha,beta LQ parameters for the actual radiation.
#' @param alphaX,betaX LQ parameters for the reference radiation.
#' @param dose dose (Gy)
#'
#' @return the value of the RBE
#'
#' @family Basic Biological Models
#' @export
rbe.evaluate <- function(alpha=alpha, beta=beta, dose=dose, alphaX=alphaX, betaX=betaX) {
R <- alphaX/betaX
RBE.alpha <- alpha/alphaX
RBE.beta2 <- beta/betaX
RBE <- R*( -1 + sqrt( 1 + 4/R * ( RBE.alpha*dose + RBE.beta2*dose^2/R ) ) ) / (2*dose)
return(RBE)
}
# GENERAL TCP/NTCP MODELS ------------------------------------------------------
#'TCP (Poissonian Model, Biological Dose)
#'
#' Evaluate the Tumor Control Probability (TCP) using the Poissonian Model from the knowledge of the biological dose distribution.
#'
#' The default parameters are taken from the example in [in Bentzen1997]
#' plus the R=alphaX/betaX parameter for tumoral tissues.
#'
#' By default it is assumed that the \code{dvh.dose} contains the biological dose dvh [Gy(E)]. Alternatively it could contain the "physical" dose [Gy] and the biological dose is evaluated by specifying the RBE per fraction (assumed to be constant over the irradiated volume).
#'
#' If alphaX and betaX are not explicity defined, they were deduced from \code{D50}, \code{gam} and \code{R}
#'
#' @param dvh.dose a dvh object in which is stored the dose (biological dose) for the specific volume.
#' @param D50,R,gam tissue specific biological parameters: tolerance dose at 50\%, R=alphaX/betaX, clonogenic density.
#' @param alphaX,betaX alternatively to (D50, R, gam), it is possible to use directly the LQ parameters as biological parameters. These parameters can be vectors.
#' @param d dose per fraction.
#' @param Nf number of fractions.
#' @param G Lea-Catcheside dose-protraction factor (if not define it is assumed complete repair between fractions, G=1/Nf).
#' @param RBE Relative Biological Effectivenes.
#'
#' @return The TCP value
#'
#' @family General TCP/NTCP Models
#' @export
#'
tcp.lqp <- function(dvh.dose=dvh.dose, D50=66, R=10, gam=2, d=2, alphaX=NULL, betaX=NULL, Nf=30, G=NULL, RBE=1)
{
# costanti
e <- exp(1)
# ricava i valori di alphaX e betaX
if((is.null(alphaX) | is.null(betaX))) {
alphaX <- alpha.beta.lqp(R, D50, gam, d)$alphaX
betaX <- alpha.beta.lqp(R, D50, gam, d)$betaX
}
#cat('LQ parameters:\n')
#cat('alphaX =', mean(alphaX), 'Gy^(-1)\n')
#cat('betaX =', mean(betaX), 'Gy^(-2)\n')
# Lea-Catcheside dose-protraction factor (assume repair completo tra frazioni)
if(is.null(G)) {G <- 1/Nf}
# volume elemetare (frazione):
dv <- 1/dvh.dose$Nvoxel
# dose assorbita per frazione
D.tot <- dvh.dose$value * Nf * RBE # dose totale assorbita
#print(length(alphaX))
S <- exp(-alphaX*D.tot - G*betaX*D.tot*D.tot)
return( exp(-exp(e*gam)*dv*sum(S)) )
}
#' NTCP (Biological Dose, Relative Seriality Model [Kallman1992a])
#'
#' It Uses as default parameters those of the "spinal cord" [from Su2010].
#'
#' @param dvh.dose a dvh object in which is stored the dose (biological dose) for the specific volume.
#' @param D50,R,gam tissue specific biological parameters: tolerance dose at 50\%, R=alphaX/betaX, clonogenic density.
#' @param alphaX,betaX alternatively to (D50, R, gam), it is possible to use directly the LQ parameters as biological parameters. These parameters can be vectors.
#' @param s relative seriality parameter of the Kallman's model.
#' @param d dose per fraction.
#' @param Nf number of fractions.
#' @param G Lea-Catcheside dose-protraction factor (if not define it is assumed complete repair between fractions, G=1/Nf).
#' @param RBE Relative Biological Effectivenes.
#'
#' @return The NTCP value
#'
#' @family General TCP/NTCP Models
#' @export
#'
ntcp.kallman <- function(dvh.dose, D50=57, R=3, gam=6.7, s=1, d=2, alphaX=NULL, betaX=NULL, Nf=30, G=NULL, RBE=1)
{
# costanti
e <- exp(1)
# ricava i valori di alphaX e betaX
if((is.null(alphaX) | is.null(betaX))) {
alphaX <- alpha.beta.lqp(R, D50, gam, d)$alphaX
betaX <- alpha.beta.lqp(R, D50, gam, d)$betaX
}
#cat('LQ parameters:\n')
#cat('alphaX =', mean(alphaX), 'Gy^(-1)\n')
#cat('betaX =', mean(betaX), 'Gy^(-2)\n')
# Lea-Catcheside dose-protraction factor (assume repair completo tra frazioni)
if(is.null(G)) {G <- 1/Nf}
# volume elemetare (frazione):
dv <- 1/dvh.dose$Nvoxel
# dose assorbita per frazione
D.tot <- dvh.dose$value * Nf * RBE # dose totale assorbita
S <- exp(-alphaX*D.tot - G*betaX*D.tot*D.tot)
return( (1 - prod( (1 - exp(-s*exp(e*gam)*S))^dv ) )^(1/s) )
}
#' TCP (Biological Dose, Poissonian Model, Variable fractions)
#'
#' TCP evaluation from the biological dose, using arbitrary set of fractions.
#'
#' The total treatment correspond to a "train" of acute doses distribution (with no time correlation, i.e. full-repair is assumed).
#' The biological dose distribution, for each voxel, can be different among fractions.
#' Hence a list of DVH objects is needed. Each element of the list correspond to a DHV evaluated for a specific fraction.
#'
#' By default it is assumed that the \code{dvh.dose} contains the biological dose dvh [Gy(E)]. Alternatively it could contain the "physical" dose [Gy] and the biological dose is evaluated by specifying the RBE per fraction (assumed to be constant over the irradiated volume).
#'
#' If alphaX and betaX are not explicity defined, they were deduced from \code{D50}, \code{gam} and \code{R}.
#'
#' @param dvhs a list of biological dose DVHs
#' @param D50,R,gam tissue specific biological parameters: tolerance dose at 50\%, R=alphaX/betaX, clonogenic density
#' @param d dose per fraction
#' @param Nf number of fractions
#' @param RBE Relative Biological Effectivenes
#'
#' @return The TCP value
#'
#' @family General TCP/NTCP Models
#' @export
#'
tcp.lqp.train <- function(dvhs, D50=66, R=10, gam=2, d=2, alphaX=NULL, betaX=NULL, RBE=1)
{
# costanti
e <- exp(1)
# ricava i valori di alphaX e betaX
if((is.null(alphaX) | is.null(betaX))) {
alphaX <- alpha.beta.lqp(R, D50, gam, d)$alphaX
betaX <- alpha.beta.lqp(R, D50, gam, d)$betaX
}
#cat('LQ parameters:\n')
#cat('alphaX =', mean(alphaX), 'Gy^(-1)\n')
#cat('betaX =', mean(betaX), 'Gy^(-2)\n')
# numero di frazioni
Nf <- length(dvhs)
#message('Number of fractions: ', Nf)
# volume elemetare (frazione). Assume che le divere frazioni abbiano la stessa griglia di voxel
dv <- 1/dvhs[[1]]$Nvoxel
# sopravvivenza finale per voxel
S <- dvhs[[1]]$value*0 + 1 # vettore di sopravvivenze finali per ogni voxel
for(i in 1:Nf) {
df <- dvhs[[i]]$value * RBE
s <- exp(-alphaX*df - betaX*df*df)
S <- S*s # assume scorrelazione temporale tra le frazioni
}
return( exp(-exp(e*gam)*dv*sum(S)) )
}
#' NTCP (Biological Dose, Relative Seriality Model [Kallman1992a], Variable fractions)
#'
#' NTCP evaluation from the biological dose, using arbitrary set of fractions.
#'
#' The total treatment correspond to a "train" of acute doses distribution (with no time correlation, i.e. full-repair is assumed).
#' The biological dose distribution, for each voxel, can be different among fractions.
#' Hence a list of DVH objects is needed. Each element of the list correspond to a DHV evaluated for a specific fraction.
#'
#' By default it is assumed that the \code{dvh.dose} contains the biological dose dvh [Gy(E)]. Alternatively it could contain the "physical" dose [Gy] and the biological dose is evaluated by specifying the RBE per fraction (assumed to be constant over the irradiated volume).
#'
#' If alphaX and betaX are not explicity defined, they were deduced from \code{D50}, \code{gam} and \code{R}.
#'
#' @param dvhs a list of biological dose DVHs
#' @param D50,R,gam tissue specific biological parameters: tolerance dose at 50\%, R=alphaX/betaX, clonogenic density
#' @param s relative seriality parameter
#' @param d dose per fraction
#' @param Nf number of fractions
#' @param RBE Relative Biological Effectivenes
#'
#' @return The NTCP value
#'
#' @family General TCP/NTCP Models
#' @export
ntcp.kallman.train <- function(dvhs, D50=57, R=3, gam=6.7, s=1, d=2, alphaX=NULL, betaX=NULL, RBE=1)
{
# costanti
e <- exp(1)
# ricava i valori di alphaX e betaX
if((is.null(alphaX) | is.null(betaX))) {
alphaX <- alpha.beta.lqp(R, D50, gam, d)$alphaX
betaX <- alpha.beta.lqp(R, D50, gam, d)$betaX
}
#cat('LQ parameters:\n')
#cat('alphaX =', mean(alphaX), 'Gy^(-1)\n')
#cat('betaX =', mean(betaX), 'Gy^(-2)\n')
# numero di frazioni
Nf <- length(dvhs)
#message('Number of fractions: ', Nf)
# volume elemetare (frazione). Assume che le divere frazioni abbiano la stessa griglia di voxel
dv <- 1/dvhs[[1]]$Nvoxel
# sopravvivenza finale per voxel
S <- dvhs[[1]]$value*0 + 1 # vettore di sopravvivenze finali per ogni voxel
for(i in 1:Nf) {
df <- dvhs[[i]]$value * RBE
sf <- exp(-alphaX*df - betaX*df*df)
S <- S*sf # assume scorrelazione temporale tra le frazioni
}
return( (1 - prod( (1 - exp(-s*exp(e*gam)*S))^dv ) )^(1/s) )
}
#' TCP (Poissonian Model)
#'
#' Evaluate TCP directly from the DVH of the "survivals" (for a given endpoint and VOI), without using the biological dose and/or indirect radiobiological parameters for a reference radiation.
#'
#' It assumes identical and uncorrelated fractions (complete-repair)
#'
#' @param dvh.survival DVH of the cell "survivals" for a single fraction
#' @param gam radiobiological parameter, related to the clonogenic density
#' @param Nf number of fractions
#'
#' @return the TCP value.
#'
#' @family General TCP/NTCP Models
#' @export
tcp.S <- function(dvh.survival=dvh.survival, gam=2, Nf=30) {
# costanti
e <- exp(1)
# volume elementare (frazione). Assume che le diverse frazioni abbiano la stessa griglia di voxel
dv <- 1/dvh.survival$Nvoxel
# Sopravvivenza complessiva (assume completa scorrelazione temporale e che le diverse frazioni siano identiche)
S <- dvh.survival$value^Nf
return( exp(-exp(e*gam)*dv*sum(S)) )
}
#' NTCP (Relative Seriality Model [Kallman1992a])
#'
#' Evaluate the NTCP directly from the DVH of the "survivals" (for a given endpoint and VOI), without using the biological dose and/or indirect radiobiological parameters for a reference radiation.
#'
#' It assumes identical and uncorrelated fractions (complete-repair)
#'
#' @param dvh.survival DVH of the cell "survivals" for a single fraction
#' @param gam radiobiological parameter, related to the clonogenic density
#' @param s relative serialiy parameter
#' @param Nf number of fractions
#'
#' @return The NTCP value.
#'
#' @family General TCP/NTCP Models
#' @export
#'
ntcp.S <- function(dvh.survival, gam=2, s=1, Nf=30) {
# costanti
e <- exp(1)
# volume elementare (frazione). Assume che le diverse frazioni abbiano la stessa griglia di voxel
dv <- 1/dvh.survival$Nvoxel
# Sopravvivenza complessiva (assume completa scorrelazione temporale e che le diverse frazioni siano identiche)
S <- dvh.survival$value^Nf
return( (1 - prod( (1 - exp(-s*exp(e*gam)*S))^dv ) )^(1/s) )
}
#LEM/MKM MODELS ----------------------------------------------------------------
#' Evaluation of alpha and beta LQ parameter (fast MKM)
#'
#' The evaluation is performed for a monoenergetic ion. It uses a C++ MKM "rapid" implementation for the evaluation. It accepts a single set of MKM parameter (\code{alphaX}, \code{betaX}, \code{rN}, \code{rd}), or alternatively a full range of variability (min,max) for each parameter.
#'
#' @param alphaX,betaX,rN,rd MKM "biological" parameter associated to a specific biological tissue/cell line:
#' \itemize{
#' \item{alphaX} LQ alpha parameter of the reference radiation [Gy^-1]
#' \item{betaX} LQ beta parameter of the reference radiation [Gy^-2]
#' \item{rN} cell nucleus radius [um]
#' \item{rd} domain radius [um]
#' }
#'
#' @param alphaX.min,alphaX.max range of variability for parameter \code{alphaX}
#' @param alphaX.N number of step for \code{alphaX}
#' @param betaX.min,betaX.max range of variability for parameter \code{betaX}
#' @param betaX.N number of step for \code{betaX}
#' @param rN.min,rN.max range of variability for parameter \code{rN}
#' @param rN.N number of step for \code{rN}
#' @param rd.min,rd.max range of variability for parameter \code{rd}
#' @param rd.N number of step for \code{rd}
#' @param cellType name of the tissue/cell line (optional)
#' @param particel type. Available ions: 'H', 'He', 'Li', 'Be', 'B,', 'C', 'N', 'O', 'F', 'Ne'.
#' @param energies vector of energies for the particle [MeV]
#' @param lets vector of LETs for the particle [keV/um]. It is used if \code{energies} is \code{NULL}.
#'
#' @return a data.frame containing all the information specified including the alpha and beta MKM evaluation (note, in the MKM implementation beta = betaX).
#'
#' @family LEM/MKM Models
#' @export
# utilizza il main main_alpha_beta_parameter_study.cc
# crea una stringa di argomenti da passare a alpha_beta_parameter_study, predisposta
# per il calcolo MKM "rapido".
#' @param calculusType the type of the evaluations options are:
#' \itemize{
#' \item{rapidKase} implementation of MKM as described in (Kase2008)
#' \item{rapidMKM} more correct implementation for the MKM from the one described in (Kase2008)
#' \item{newMKMKase} "rapidKase" implementation + non poissionian correction for beta
#' \item{newMKM} "rapidMKM" implementation + non poissionian correction for beta
#' \item{slow_alphaIon_betaIon} slow Monte Carlo evaluation (compatible with MKM, LEMI, LEMII and LEMIII).
#' }
# In output ritorna un data.frame con i parametri corrispondenti. I valori di
# output vengono passati attraverso il file temporaneo salvato dall'eseguibile, che
# viene cancellato subito dopo. Il file temporaneo è identificato dal nome della
# cellType. Nota: potrebbe essere possibile creare automaticamente N cellType
# e quindi lanciare in parallelo N processi di calcolo...?
alpha.beta.mkm <- function(alphaX=0.1295, betaX=0.03085, rN=4, rd=0.31,
alphaX.min=NULL, alphaX.max=NULL, alphaX.N=NULL,
betaX.min=NULL, betaX.max=NULL, betaX.N=NULL,
rN.min=NULL, rN.max=NULL, rN.N=NULL,
rd.min=NULL, rd.max=NULL, rd.N=NULL,
cellType=NULL,
particleType='H',
energies=NULL, lets=NULL,
calculusType='rapidMKM',
ignore.stdout=TRUE, ignore.stderr=TRUE)
{
model <- 'MKM'
precision <- 1
#lem.setenv=get('lem.setenv', envir=dektoolsEnv)
# nome cellType
if(is.null(cellType)) {cellType <- paste('R', sprintf("%06d", round(runif(1,min=0,max=1e6))), sep='')}
# valori dei parametri MKM
if (!is.null(alphaX)) {a <- paste(alphaX, alphaX, 1)} else {a <- paste(alphaX.min, alphaX.max, alphaX.N)}
if (!is.null(betaX)) {b <- paste(betaX, betaX, 1)} else {b <- paste(betaX.min, betaX.max, betaX.N)}
if (!is.null(rN)) {n <- paste(rN, rN, 1)} else {n <- paste(rN.min, rN.max, rN.N)}
if (!is.null(rd)) {d <- paste(rd, rd, 1)} else {d <- paste(rd.min, rd.max, rd.N)}
# check per vedere se il calcolo è sulle energie o sul let (e controllo estremi):
if(length(energies)>0) {
energyType <- 'energy'
e <- paste(energies, collapse=' ')
} else if(length(lets)>0) {
lets <- lets[lets >= srim.let.min[particleType] & lets <= srim.let.max[particleType]]
energyType <- 'let'
e <- paste(lets, collapse=' ')
} else if(length(energies)==0 & length(lets)==0){
stop('energies/lets not specified.')
}
# output file
of <- paste(particleType, cellType, model, calculusType, sep='_')
of <- paste(of, 'csv', sep='.')
# costruzione linea di comando:
s.args <- paste(cellType, model, calculusType, a, b, n, d, particleType, precision, energyType, e)
# cmd <- paste('.', lem.setenv, '; survival_alpha_beta_parameter_study', s.args)
cmd <- paste('survival', s.args)
if(!ignore.stdout) {message(cmd)}
t <- system.time(system(cmd, ignore.stdout=ignore.stdout, ignore.stderr=ignore.stderr))
#message('time elapsed: ', t)
# legge file temporaneo salvato
out.df <- read.csv(of)
# cancella file temporaneo
system(paste('rm', of))
return(out.df)
}
#' Evaluation of ion alpha and beta LQ parameters
#'
#' The evaluation is performed for a monoenergetic ion. It uses a C++ implementation for the evaluation. It accepts a single set of parameter depending on the choosen radiobiological model, or alternatively a full range of variability (min,max) for each parameter.
#'
#' @param model the name of the model (options: 'MKM', 'LEMI', 'LEMII', 'LEMIII')
#' @param model.parameters a dataframe containing the "biological" parameter associated to a specific biological tissue/cell line. In the case of the MKM:
#' \itemize{
#' \item{alpha0} LQ alpha parameter of the reference radiation [Gy^-1]
#' \item{beta0} LQ beta parameter of the reference radiation [Gy^-2]
#' \item{rN} cell nucleus radius [um]
#' \item{rd} domain radius [um]
#' }
#' In the case of LEMI, LEMII, LEMIII:
#' \itemize{
#' \item{alpha0} LQ alpha parameter of the reference radiation [Gy^-1]
#' \item{beta0} LQ beta parameter of the reference radiation [Gy^-2]
#' \item{rN} cell nucleus radius [um]
#' \item{Dt} threshold dose [Gy]
#' }
#' Optionally, instead of the specific values, the range of the parameters can be specified in the same data.frame. In the case of MKM:
#' \itemize{
#' \item{alpha0.min, alpha0.max, alpha0.N} range of variability for parameter \code{alpha0} and number of steps
#' \item{beta0.min, beta0.max, beta0.N} range of variability for parameter \code{beta0} and number of steps
#' \item{rN.min, rN.max, rN.N} range of variability for parameter \code{rN} and number of steps
#' \item{rd.min, rN.max} domain radius [um]
#' }
#'
#' @param cellType name of the tissue/cell line (optional)
#' @param particle available particles: 'H', 'He', 'Li', 'Be', 'B,', 'C', 'N', 'O', 'F', 'Ne'.
#' @param energies vector of energies for the particle [MeV]
#' @param lets vector of LETs for the particle [keV/um]. It is used if \code{energies} is \code{NULL}.
#' @param calculusType the type of the evaluations options are:
#' \itemize{
#' \item{rapidKase} implementation of MKM as described in (Kase2008)
#' \item{rapidMKM} more correct implementation for the MKM from the one described in (Kase2008)
#' \item{newMKMKase} "rapidKase" implementation + non poissionian correction for beta
#' \item{newMKM} "rapidMKM" implementation + non poissionian correction for beta
#' \item{rapidScholz} fast implementation of the LEMI, LEMII, LEMIII (Scholz2000)
#' \item{rapidRusso} fast implementation of the LEMI, LEMII, LEMIII, more accurate than the \code{rapidScholz} (Russo2010)
#' \item{slow_alphaIon_betaIon} slow Monte Carlo evaluation (compatible with MKM, LEMI, LEMII and LEMIII).
#' }
#' @param precision used only for the Monte Carlo evaluation. If precision < 1.0 it represents the target relative standard deviation to be reached in the evaluations of the cell survivals. If precision >= 1 it represents the number of cells to be simulated.
#' @param get.raw.data If calculusType=slow_alphaIon_betaIon, returns the raw data (survival vs dose for each simulated cell) from which alpha and beta parameters can be extracted, otherwhise it is ignored.
#'
#' @return a data.frame containing all the information specified including the alpha and beta MKM evaluation (note, in the MKM implementation beta = betaX).
#'
#' @family LEM/MKM Models
#' @export
alpha.beta.ion.range <- function(model='MKM',
model.parameters=data.frame(alpha0=0.1295, beta0=0.03085, rN=4, rd=0.31),
cell.name=NULL,
particle='H',
energies=NULL, lets=NULL,
calculusType='rapidMKM', precision=0.5,
ignore.stdout=TRUE, ignore.stderr=TRUE,
get.raw.data=FALSE, remove.temp.files=TRUE)
{
# nome cellType
if(is.null(cell.name)) {cell.name <- paste('R', sprintf("%06d", round(runif(1,min=0,max=1e6))), sep='')}
# nomi dei parametri
if(model=='MKM') {
parameter.names <- c('alpha0', 'beta0', 'rN', 'rd')
parameter.names.min <- c('alpha0.min', 'beta0.min', 'rN.min', 'rd.min')
parameter.names.max <- c('alpha0.max', 'beta0.max', 'rN.max', 'rd.max')
parameter.names.N <- c('alpha0.N', 'beta0.N', 'rN.N', 'rd.N')
}
else if(model=='LEMI' | model=='LEMII' | model=='LEMIII') {
parameter.names <- c('alpha0', 'beta0', 'rN', 'Dt')
parameter.names.min <- c('alpha0.min', 'beta0.min', 'rN.min', 'Dt.min')
parameter.names.max <- c('alpha0.max', 'beta0.max', 'rN.max', 'Dt.max')
parameter.names.N <- c('alpha0.N', 'beta0.N', 'rN.N', 'Dt.N')
} else {
stop('Model', model, 'not defined.')
}
# recupera parametri
if(all(parameter.names %in% names(model.parameters))) {
p <- model.parameters[parameter.names]
p.min <- p.max <- p.N <- NULL
} else if(all(c(parameter.names.min, parameter.names.max, parameter.names.N) %in% names(model.parameters))) {
p.min <- model.parameters[parameter.names.min];
p.max <- model.parameters[parameter.names.max];
p.N <- model.parameters[parameter.names.N];
p <- NULL
} else {
stop('Parameter not well defined. (', names(model.parameters), ')')
}
# valori dei parametri da passare a pure-survival
if (!is.null(p[1])) {a <- paste(p[1], p[1], 1)} else {a <- paste(p.min[1], p.max[1], p.N[1])}
if (!is.null(p[2])) {b <- paste(p[2], p[2], 1)} else {b <- paste(p.min[2], p.max[2], p.N[2])}
if (!is.null(p[3])) {n <- paste(p[3], p[3], 1)} else {n <- paste(p.min[3], p.max[3], p.N[3])}
if (!is.null(p[4])) {d <- paste(p[4], p[4], 1)} else {d <- paste(p.min[4], p.max[4], p.N[4])}
# check per vedere se il calcolo è sulle energie o sul let (e controllo estremi):
if(length(energies)>0) {
energyType <- 'energy'
e <- paste(energies, collapse=' ')
} else if(length(lets)>0) {
lets <- lets[lets >= srim.let.min[particle] & lets <= srim.let.max[particle]]
energyType <- 'let'
e <- paste(lets, collapse=' ')
} else if(length(energies)==0 & length(lets)==0){
stop('energies/lets not specified.')
}
# output file
of <- paste(particle, cell.name, model, calculusType, sep='_')
of <- paste(of, 'csv', sep='.')
# creazione cartella per calcolo MC
if(calculusType=='slow_alphaIon_betaIon') {dir.create(paste0('survival-data-', particle), showWarnings = FALSE)}
# costruzione linea di comando:
s.args <- paste(cell.name, model, calculusType, a, b, n, d, particle, precision, energyType, e)
# cmd <- paste('.', lem.setenv, '; survival_alpha_beta_parameter_study', s.args)
cmd <- paste('survival', s.args)
if(!ignore.stdout) {message(cmd)}
t <- system.time(system(cmd, ignore.stdout=ignore.stdout, ignore.stderr=ignore.stderr))
#message('time elapsed: ', t)
# legge file temporaneo salvato
out.df <- read.csv(of)
out.df$cell <- cell.name
out.df$particle <- particle
# cancella file temporaneo
if(remove.temp.files) {
#system(paste('rm', of))
unlink(of)
if(calculusType=='slow_alphaIon_betaIon') {
unlink(paste0('survival-data-', particle), recursive = TRUE)
}
}
return(out.df)
}
#' Evaluation of ion alpha and beta LQ parameters
#'
#' The evaluation is performed for a monoenergetic ion. It uses a C++ implementation for the evaluation. It accepts a sequence of parameters and ion specifications stored in an input data.frame.
#'
#' @param model the name of the model (options: 'MKM', 'LEMI', 'LEMII', 'LEMIII')
#' @param calculusType the type of the evaluations options are:
#' \itemize{
#' \item{rapidKase} implementation of MKM as described in (Kase2008)
#' \item{rapidMKM} more correct implementation for the MKM from the one described in (Kase2008)
#' \item{newMKMKase} "rapidKase" implementation + non poissionian correction for beta
#' \item{newMKM} "rapidMKM" implementation + non poissionian correction for beta
#' \item{rapidScholz} fast implementation of the LEMI, LEMII, LEMIII (Scholz2000)
#' \item{rapidRusso} fast implementation of the LEMI, LEMII, LEMIII, more accurate than the \code{rapidScholz} (Russo2010)
#' \item{slow_alphaIon_betaIon} slow Monte Carlo evaluation (compatible with MKM, LEMI, LEMII and LEMIII).
#' }
#' @param parameters input data.frame containing the cell name (column 'cell', optional) "biological" parameters associated to a specific biological tissue/model, and the specification of the primary particles (ion type and energy and/or LET). All the parameters have to be put in a dataframe. Each row of the dataframe represent the irradiation of a specific tissue/model with a specific particle. In the case of the MKM the biological parameters
#' \itemize{
#' \item{alpha0} LQ alpha parameter of the reference radiation [Gy^-1]
#' \item{beta0} LQ beta parameter of the reference radiation [Gy^-2]
#' \item{rN} cell nucleus radius [um]
#' \item{rd} domain radius [um]
#' }
#' In the case of LEMI, LEMII, LEMIII:
#' \itemize{
#' \item{alpha0} LQ alpha parameter of the reference radiation [Gy^-1]
#' \item{beta0} LQ beta parameter of the reference radiation [Gy^-2]
#' \item{rN} cell nucleus radius [um]
#' \item{Dt} threshold dose [Gy]
#' }
#' The ion is specified with the columns:
#' \itemize{
#' \item{particle} ion type: 'H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F' or 'Ne'
#' \item{energy} Energy of the particle in MeV/u
#' \item{let} LET of the particle in keV/um. It is used only if the energy is not specified.
#' }
#' #' @family LEM/MKM Models
#' @export
#' @import data.table
alpha.beta.ion.sequence <- function(model='MKM',
parameters=data.frame(cell='test', alpha0=0.1295, beta0=0.03085, rN=4, rd=0.31,
particle='H', energy=50, let=NULL),
calculusType='rapidMKM', precision=0.5,
ignore.stdout=TRUE, ignore.stderr=TRUE,
remove.temp.files=TRUE)
{
#parameters <- data.table(parameters)
parameters.column.index <- which(names(parameters) %in% c('alpha0', 'beta0', 'rN', 'rd', 'Dt'))
if(length(parameters.column.index)!=4) {stop('alpha.beta.ion.sequence: error in the definition of the parameters')}
# dataframe di output
alpha.beta.out <- NULL
# loop su id.par
Npar <- nrow(parameters)
for(id in 1:Npar) {
if(nrow(parameters)==0) {break}
par <- parameters[1, parameters.column.index]
index.par <- parameters[,parameters.column.index[1]] == par[1,1] &
parameters[,parameters.column.index[2]] == par[1,2] &
parameters[,parameters.column.index[3]] == par[1,3] &
parameters[,parameters.column.index[4]] == par[1,4]
parameters.m <- parameters[index.par, ]
cell <- parameters.m$cell[1]
message('cell: ', cell)
print(par)
#print(summary(parameters.m))
# loop su particles
particles <- unique(parameters.m$particle)
for(ip in 1:length(particles)) {
message('particle: ', particles[ip])
parameters.p <- subset(parameters.m, particle==particles[ip])
energies <- parameters.p$energy
lets <- parameters.p$let
alpha.beta.tmp <- alpha.beta.ion.range(model=model, model.parameters=par, cell.name=cell, particle=particles[ip], energies=energies, lets=lets, calculusType=calculusType, precision=precision, ignore.stdout=ignore.stdout, ignore.stderr=ignore.stderr, remove.temp.files=remove.temp.files)
alpha.beta.out <- rbind(alpha.beta.out, alpha.beta.tmp)
}
parameters <- parameters[-which(index.par),]
}
return(alpha.beta.out)
}
#' @export
#' @import data.table
alpha.beta.ion.sequence2 <- function(model='MKM',
parameters=data.frame(cell='test', alpha0=0.1295, beta0=0.03085, rN=4, rd=0.31,
particle='H', energy=50, let=NA),
calculusType='rapidMKM', precision=0.5,
ignore.stdout=TRUE, ignore.stderr=TRUE,
remove.temp.files=TRUE)
{
#parameters <- data.table(parameters)
parameters.column.index <- which(names(parameters) %in% c('alpha0', 'beta0', 'rN', 'rd', 'Dt'))
if(length(parameters.column.index)!=4) {stop('alpha.beta.ion.sequence: error in the definition of the parameters')}
# fa diventare data.table....
#if('data.table' %in% class(parameters)) {
# message('casting in data.table')
# parameters <- as.data.table(parameters)
#}
# aggiunge identificativo unico per i parametri radiobiologici
#parameters$id.par <- interaction(parameters[, parameters.column.index])
# ho perso le chiavi
#setkey(parameters, id.par)
# dataframe di output
alpha.beta.out <- NULL
# loop su id.par
message('identifying models...')
id.pars <- unique(parameters[, parameters.column.index])
message(nrow(id.pars), ' different models')
for(id in 1:nrow(id.pars)) {
parameters.m <- parameters[parameters[,parameters.column.index[1]] == id.pars[id,1] &
parameters[,parameters.column.index[2]] == id.pars[id,2] &
parameters[,parameters.column.index[3]] == id.pars[id,3] &
parameters[,parameters.column.index[4]] == id.pars[id,4], ]
biological.parameters <- id.pars[id,]
cell <- parameters.m$cell[1]
message('cell: ', cell)
print(biological.parameters)
# loop su particles
particles <- unique(parameters.m$particle)
for(ip in 1:length(particles)) {
message('particle: ', particles[ip])
parameters.p <- subset(parameters.m, particle==particles[ip])
energies <- parameters.p$energy
lets <- parameters.p$let
alpha.beta.tmp <- alpha.beta.ion.range(model=model, model.parameters=biological.parameters, cell.name=cell, particle=particles[ip], energies=energies, lets=lets, calculusType=calculusType, precision=precision, ignore.stdout=ignore.stdout, ignore.stderr=ignore.stderr, remove.temp.files=remove.temp.files)
alpha.beta.out <- rbind(alpha.beta.out, alpha.beta.tmp)
}
}
return(alpha.beta.out)
}
#' Evaluate alpha function (MKM)
#'
#' Returns a function alpha(LET) or alpha(specific energy), using the MKM model. The functions is an interpolation functions that uses a set of alpha values evaluated over the defined set of LETs (or specific energies).
#'
#' @param alphaX,betaX,rN,rd The radiobiological parameter of the MKM.
#' @param cellType Optionally you can associate a cell name tag (I have to figure yet what use to do with that).
#' @param particleType The particle specification ('H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F')
#' @param energies,lets A vector of specific energies (MeV/u) or LETs (keV/um).
#' @param use.limits Filters the energies or the LETs values to rule out values too high or too low.
#' @param calculusType the type of the evaluations options are:
#' \itemize{
#' \item{rapidKase} implementation of MKM as described in (Kase2008)
#' \item{rapidMKM} more correct implementation for the MKM from the one described in (Kase2008)
#' \item{newMKMKase} "rapidKase" implementation + non poissionian correction for beta
#' \item{newMKM} "rapidMKM" implementation + non poissionian correction for beta
#' \item{slow_alphaIon_betaIon} slow Monte Carlo evaluation.
#' }
#'
#' @family LEM/MKM Models
#' @export
#'
alpha.fun.mkm <- function(alphaX=0.1295, betaX=0.03085, rN=4, rd=0.31,
cellType=NULL,
particleType='H',
energies=NULL, lets=NULL,
use.limits=FALSE, outmessages=FALSE,
calculusType='rapidMKM')
{
model <- 'MKM'
calculusType <- 'rapidMKM'
precision <- 1
#lem.setenv=get('lem.setenv', envir=dektoolsEnv)
# nome cellType
if(is.null(cellType)) {
cellType <- paste('R', sprintf("%06d", round(runif(1,min=0,max=1e6))), sep='')
#message('cellType = ', cellType)
}
# valori dei parametri MKM
a <- paste(alphaX, alphaX, 1)
b <- paste(betaX, betaX, 1)
n <- paste(rN, rN, 1)
d <- paste(rd, rd, 1)
# check per vedere se il calcolo è sulle energie o sul let:
if(!is.null(energies)) {
energyType <- 'energy'
e <- paste(energies, collapse=' ')
ee <- energies
} else {
lets <- lets[lets >= srim.let.min[particleType] & lets <= srim.let.max[particleType]]
energyType <- 'let';
e <- paste(lets, collapse=' ')
ee <- lets
}
# check per vedere se siamo fuori limiti (restituisce valori alti di alpha, da usare nelle ottimizzazioni)
if(use.limits) {
#message('using limits')
if(!prod(c(alphaX, betaX, rN, rd)>mkm.lower & c(alphaX, betaX, rN, rd)<mkm.upper)) {
mkm.fun <- approxfun(ee, rep(1e8, length(ee)), method='linear')
#message(alphaX, betaX, rN, rd, ' out of limits!', sep=' ')
return(mkm.fun)
}
}
# output file
of <- paste(particleType, cellType, model, calculusType, sep='_')
of <- paste(of, 'csv', sep='.')
# costruzione linea di comando:
s.args <- paste(cellType, model, calculusType, a, b, n, d, particleType, precision, energyType, e)
#cmd <- paste('.', lem.setenv, '; survival_alpha_beta_parameter_study', s.args)
cmd <- paste('survival', s.args)
#print(cmd)
if(outmessages) {ignore.stdout=FALSE; ignore.stderr=FALSE} else {ignore.stdout=TRUE; ignore.stderr=TRUE}
t <- system.time(system(cmd, ignore.stdout=TRUE, ignore.stderr=TRUE))
#message('time elapsed: ', t)
# legge file temporaneo salvato
out.df <- read.csv(of)
# cancella file temporaneo
system(paste('rm', of))
# crea funzione interpolazione
mkm.fun <- approxfun(ee, out.df$alpha, method='linear')
#print(c(alphaX, betaX, rN, rd))
return(mkm.fun)
}
# VALUES BIOLOGICAL EVALUATIONS ------------------------------------------------
#' Calcola la sopravvivenza secondo il modello LM e la aggiunge in values
#'
#' m, c e beta sono i parametri del modello
#' m -> coefficiente angolare
#' c -> costante
#' beta -> valore di beta (indipendente dal LETd)
#' Nota: l'oggetto values deve avere le due variabili 'Dose[Gy]' e 'DoseAveragedLET[keV/um]'
#'
#' @family Values Bioogical Evaluations
#' @export
#'
values.add.survival.LM <- function(values, alpha0, m, beta, add.alpha=TRUE)
{
# identifica gli indici per 'Dose[Gy]' e 'DoseAveragedLET[keV/um]'
index.dose<- which(values$variables=='Dose[Gy]')
index.LETd <- which(values$variables=='DoseAveragedLET[keV/um]')
# estrae dose e LETd da values
dose <- values$values[index.dose,,,]
LETd <- values$values[index.LETd,,,]
# calcola un nuovo array 'Survival.LM' da 'Dose[Gy]' e 'DoseAveragedLET[keV/um]' usando come parametri m, c, e beta
alpha <- m*LETd + alpha0
new.array <- exp(-alpha*dose - beta*dose^2)
# aggiungi l'array calcolato all'oggetto values
values <- add.array.values(values, new.array, variable='Survival.LM')
# aggiungi opzionalmente anche l'array alpha
if(add.alpha) {values <- add.array.values(values, alpha, variable='Alpha.LM[Gy^(-1)]')}
# ritorna values
return(values)
}
#' Calcola l'RBE.alpha e l'aggiunge in values
#'
#' viene passato come argomento il modello:
#' (NULL -> usa "Alpha[Gy^(-1)]")
#' (LM -> usa "Alpha.LM[Gy^(-1)]")
#' (cMKM -> usa "Alpha.cMKM[Gy(-1)]")
#'
#' assume che la variabile alpha sia già presente in values, atrimenti la calcola
#' esplicitamente.
#'
#' nota: nell'oggetto values NON deve essere già presente la variabile "Survival"
#'
#' nota: nei parametri del modello deve essere esplcitata obbligatoriamente anche il beta
#' anche se non viene esplicitamente usatato, viene usato per calcolare le sopravvivenze
#'
#' @family Values Bioogical Evaluations
#' @export
#'
values.add.RBE.alpha <- function(values, model=NULL, alphaX=0.2, model.LM=NULL, model.cMKM=NULL)
{
if(is.null(model)) {var.alpha <- 'Alpha[Gy^(-1)]'; var.rbe <- 'RBE.alpha'}
else if(model=='LM') {var.alpha <- 'Alpha.LM[Gy^(-1)]'; var.rbe <- 'RBE.alpha.LM'}
else if(model=='cMKM') {var.alpha <- 'Alpha.cMKM[Gy^(-1)]'; var.rbe <- 'RBE.alpha.cMKM'}
v.alpha <- which(values$variables==var.alpha)
# calcola alpha ex novo
if(length(v.alpha)==0) {
if(model=='LM') {
message('Evaluating RBE.alpha using Model: LM')
values <- values.add.survival.LM(values, alpha0=model.LM$alpha0, m=model.LM$m, beta=model.LM$beta)
} else if(model=='cMKM') {
message('Evaluating RBE.alpha using Model: LM')
values <- values.add.survivals.LM(values,
alphaX=model.cMKM$alphaX, betaX=model.cMKM$betaX,
rN=model.cMKM$rN, rd=model.cMKM$rd,
beta=model.cMKM$beta,
particleType=as.charcater(model.cMKM$particleType),
beta=beta, add.alpha=TRUE)
}
}
v.alpha <- which(values$variables==var.alpha)
rbe.array <- values$values[v.alpha,,,]/alphaX
values <- add.array.values(values=values, new.array=rbe.array, variable=var.rbe)
return(values)
}
#' Calcola la sopravvivenza secondo il modello cMKM e la aggiunge in values
#'
#' alphaX, betaX, rN, rd sono i parametri del modello MKM
#' particleType -> ione di riferimento
#' beta -> valore di beta (indipendente dal LETd)
#' Nota: l'oggetto values deve avere le due variabili 'Dose[Gy]' e 'DoseAveragedLET[keV/um]'
#'
#' @family Values Bioogical Evaluations
#' @export
#'
values.add.survival.cMKM <- function(values=values, alphaX=alphaX, betaX=betaX, rN=rN, rd=rd, particleType=particleType, beta=beta, add.alpha=TRUE)
{
# identifica gli indici per 'Dose[Gy]' e 'DoseAveragedLET[keV/um]'
index.dose<- which(values$variables=='Dose[Gy]')
index.LETd <- which(values$variables=='DoseAveragedLET[keV/um]')
# estrae dose e LETd da values
dose <- values$values[index.dose,,,]
LETd <- values$values[index.LETd,,,]
# vettore di lets da interpolare
lets <- seq(min(LETd, na.rm=TRUE), max(LETd, na.rm=TRUE), length.out=30)
# calcola un nuovo array 'Survival.cMKM' da 'Dose[Gy]' e 'DoseAveragedLET[keV/um]'
alpha.mkm <- alpha.fun.mkm(alphaX=alphaX, betaX=betaX, rN=rN, rd=rd,
particleType=particleType,
lets=lets,
use.limits=FALSE)
alpha <- alpha.mkm(LETd)
new.array <- exp(-alpha*dose - beta*dose^2)
# aggiungi l'array calcolato all'oggetto values
values <- add.array.values(values=values, new.array=new.array, variable='Survival.cMKM')
# aggiungi opzionalmente anche l'array alpha
if(add.alpha) {values <- add.array.values(values, alpha, variable='Alpha.cMKM[Gy^(-1)]')}
# ritorna values
return(values)
}
#' calcola l'RBE e l'agiunge in values a partire da alpha e beta e dose
#'
#' @family Values Bioogical Evaluations
#' @export
#'
values.add.RBE <- function(values, model=NULL, alphaX=0.2, betaX=0.02, model.LM=NULL, model.cMKM=NULL, biological.dose=FALSE)
{
if(is.null(model)) {var.alpha <- 'Alpha[Gy^(-1)]'; var.beta <- 'Beta[Gy^(-2)]'; var.rbe <- 'RBE'}
else if(model=='LM') {var.alpha <- 'Alpha.LM[Gy^(-1)]'; var.beta <- 'Beta.LM[Gy^(-2)]'; var.rbe <- 'RBE.LM'}
else if(model=='cMKM') {var.alpha <- 'Alpha.cMKM[Gy^(-1)]'; var.beta <- 'Beta.cMKM[Gy^(-2)]'; var.rbe <- 'RBE.cMKM'}
v.alpha <- which(values$variables==var.alpha)
v.beta <- which(values$variables==var.beta)
# calcola alpha e beta ex novo, e quindi l'rbe
if(length(v.alpha)+length(v.beta)==0) {
# modello LM
if(model=='LM') {
message('Evaluating RBE using Model: LM')
index.dose <- which(values$variables=='Dose[Gy]')
index.LETd <- which(values$variables=='DoseAveragedLET[keV/um]')
LETd <- values$values[index.LETd,,,]
dose <- values$values[index.dose,,,]
alpha <- model.LM$m*LETd + model.LM$alpha0
beta <- alpha*0 + model.LM$beta
rbe <- rbe.evaluate(alpha=alpha, beta=beta, dose=dose, alphaX=alphaX, betaX=betaX)
values <- add.array.values(values=values, new.array=rbe, variable='RBE.LM')
if(biological.dose) {
dbio <- dose*rbe
values <- add.array.values(values=values, new.array=dbio, variable='BiologicalDose.LM[Gy(RBE)]')
}
}
# modello cMKM
else if(model=='cMKM') {
message('Evaluating RBE using Model: cMKM')
index.dose <- which(values$variables=='Dose[Gy]')
index.LETd <- which(values$variables=='DoseAveragedLET[keV/um]')
LETd <- values$values[index.LETd,,,]
dose <- values$values[index.dose,,,]
lets <- seq(min(LETd, na.rm=TRUE), max(LETd, na.rm=TRUE), length.out=30)
alpha.mkm <- alpha.fun.mkm(alphaX=model.cMKM$alphaX, betaX=model.cMKM$betaX, rN=model.cMKM$rN, rd=model.cMKM$rd,
cellType='RCell', particleType=as.character(model.cMKM$particleType), lets=lets, use.limits=FALSE)
alpha <- alpha.mkm(LETd)
beta <- alpha*0 + model.cMKM$beta
rbe <- rbe.evaluate(alpha=alpha, beta=beta, dose=dose, alphaX=alphaX, betaX=betaX)
values <- add.array.values(values=values, new.array=rbe, variable='RBE.cMKM')
if(biological.dose) {
dbio <- dose*rbe
values <- add.array.values(values=values, new.array=dbio, variable='BiologicalDose.cMKM[Gy(RBE)]')
}
}
}
return(values)
}
# LM/cMKM TCP/NTCP MODELS ------------------------------------------------------
#' calcola TCP con modello LM
#'
#' @family LM/cMKM TCP/NTCP Models
#' @export
#'
TCP.LM <- function(values, vois, voi='PTV', alpha0, beta, m, G=G, Nf=Nf){
values.S <- values.add.survival.LM(values=values, alpha0=alpha0, m=m, beta=beta)
dvh.surv <- dvh.evaluate(values=values.S, variable='Survival.LM', voi=voi, vois=vois)
TCP <- tcp.S(dvh.survival=dvh.surv, gam=G, Nf=Nf)
return(TCP)
}
#' calcola TCP con modello cMKM
#' @family LM/cMKM TCP/NTCP Models
#' @export
#'
TCP.cMKM <- function(values=values, vois=vois, voi='PTV', alphaX=alphaX, betaX=betaX, rN=rN, rd=rd, beta=beta, particleType=particleType, G=G, Nf=Nf){
values.S <- values.add.survival.cMKM(values=values, alphaX=alphaX, betaX=betaX, rN=rN, rd=rd, beta=beta, particleType=particleType)
dvh.surv <- dvh.evaluate(values=values.S, variable='Survival.cMKM', voi=voi, vois=vois)
TCP <- tcp.S(dvh.survival=dvh.surv, gam=G, Nf=Nf)
return(TCP)
}
#' calcola NTCP con modello LM
#' @family LM/cMKM TCP/NTCP Models
#' @export
#'
NTCP.LM <- function(values, vois, voi='PTV', alpha0, beta, m, G=G, s=s, Nf=Nf){
values.S <- values.add.survival.LM(values=values, alpha0=alpha0, m=m, beta=beta)
dvh.surv <- dvh.evaluate(values=values.S, variable='Survival.LM', voi=voi, vois=vois)
NTCP <- ntcp.S(dvh.survival=dvh.surv, gam=G, s=s, Nf=Nf)
return(NTCP)
}
#' calcola NTCP con modello cMKM
#'
#' @family LM/cMKM TCP/NTCP Models
#' @export
#'
NTCP.cMKM <- function(values=values, vois=vois, voi='PTV', alphaX=alphaX, betaX=betaX, rN=rN, rd=rd, beta=beta, particleType=particleType, G=G, s=s, Nf=Nf){
values.S <- values.add.survival.cMKM(values=values, alphaX=alphaX, betaX=betaX, rN=rN, rd=rd, beta=beta, particleType=particleType)
dvh.surv <- dvh.evaluate(values=values.S, variable='Survival.cMKM', voi=voi, vois=vois)
NTCP <- ntcp.S(dvh.survival=dvh.surv, gam=G, s=s, Nf=Nf)
return(NTCP)
}
# LIKELIHOOD EVALUATIONS -------------------------------------------------------
#' calcola la log likelihood di cura (TCP), per un singolo caso, rispetto al dato clinico
#'
#' Utilizza il modello lineare (alpha=alpha0 + m*LETd, beta=cost.).
#' Assume frazioni indipendenti e costanti.
#' Assume la presenza di matrice di dose fisica e di LETd per il calcolo LM.
#'
#' Np -> numero di pazienti
#' R -> numero di pazienti responder (R=Np*TCP)
#' Nf -> numero di frazioni
#'
#' @family Likelihood Evaluations
#' @export
#'
TCP.LL.LM <- function(values, vois, voi='PTV', alpha0, beta, m, G=G, Np=Np, R=R, Nf=Nf){
values.S <- values.add.survival.LM(values=values, alpha0=alpha0, m=m, beta=beta)
dvh.surv <- dvh.evaluate(values=values.S, variable='Survival.LM', voi=voi, vois=vois)
TCP <- tcp.S(dvh.survival=dvh.surv, gam=G, Nf=Nf)
LL <- log(TCP)*R+log(1-TCP)*(Np-R)
return(LL)
}
#' calcola la log likelihood di cura (TCP), per un singolo caso, rispetto al dato clinico
#'
#' Utilizza il modello cMKM.
#' Assume frazioni indipendenti e costanti.
#' Assume la presenza di matrice di dose fisica e di LETd per il calcolo LM.
#'
#' Np -> numero di pazienti
#' R -> numero di pazienti responder (R=Np*TCP)
#' Nf -> numero di frazioni
#' @family Likelihood Evaluations
#' @export
#'
TCP.LL.cMKM <- function(values, vois, voi='PTV',
alphaX, betaX, rN, rd, beta, particleType,
G, Np, R, Nf){
values.S <- values.add.survival.cMKM(values=values, alphaX=alphaX, betaX=betaX, rN=rN, rd=rd, particleType=particleType, beta=beta)
dvh.surv <- dvh.evaluate(values=values.S, variable='Survival.LM', voi=voi, vois=vois)
TCP <- tcp.S(dvh.survival=dvh.surv, gam=G, Nf=Nf)
LL <- log(TCP)*R+log(1-TCP)*(Np-R)
return(LL)
}
#' calcola la log likelihood di complicazioni (NTCP), per un singolo caso, rispetto al dato clinico
#'
#' Utilizza il modello lineare (alpha=alpha0 + m*LETd, beta=cost.).
#' Assume frazioni indipendenti e costanti.
#' Assume la presenza di matrice di dose fisica e di LETd per il calcolo LM.
#'
#' Np -> numero di pazienti
#' R -> numero di pazienti responder (R=Np*NTCP)
#' Nf -> numero di frazioni
#' @family Likelihood Evaluations
#' @export
#'
NTCP.LL.LM <- function(values, vois, alpha0, beta, m, G=G, s=s, Np=Np, R=R, Nf=Nf, voi=voi){
values.S <- values.add.survival.LM(values=values, alpha0=alpha0, m=m, beta=beta)
dvh.surv <- dvh.evaluate(values=values.S, variable='Survival.LM', voi=voi, vois=vois)
NTCP <- ntcp.S(dvh.survival=dvh.surv, gam=G, s=s, Nf=Nf)
LL <- log(NTCP)*R+log(1-NTCP)*(Np-R)
if(is.infinite(LL)) {message('NTCP: ', NTCP)}
return(LL)
}
#' calcola la log likelihood di complicazioni (NTCP), per un singolo caso, rispetto al dato clinico
#'
#' Utilizza il modello cMKM.
#' Assume frazioni indipendenti e costanti.
#' Assume la presenza di matrice di dose fisica e di LETd per il calcolo LM.
#'
#' Np -> numero di pazienti
#' R -> numero di pazienti responder (R=Np*NTCP)
#' Nf -> numero di frazioni
#' @family Likelihood Evaluations
#' @export
#'
NTCP.LL.cMKM <- function(values, vois,
alphaX, betaX, rN, rd, beta, particleType,
G=G, s=s, Np=Np, R=R, Nf=Nf, voi=voi){
values.S <- values.add.survival.cMKM(values=values, alphaX=alphaX, betaX=betaX, rN=rN, rd=rd, particleType=particleType, beta=beta)
dvh.surv <- dvh.evaluate(values=values.S, variable='Survival.LM', voi=voi, vois=vois)
NTCP <- ntcp.S(dvh.survival=dvh.surv, gam=G, s=s, Nf=Nf)
LL <- log(NTCP)*R+log(1-NTCP)*(Np-R)
if(is.infinite(LL)) {message('NTCP: ', NTCP)}
return(LL)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.