R/SPRModel_Langmuir.R

In ffeng23/ADASPR: An R package to characterize the distributed affinities in a multi-ligand system

#specify the dependency
#' @include SPRData.R
NULL

				
#'@title S4 class of the SPR Langmuir model
#'
#'@description Langmuir model object
#'
#'@details check the manual for the detailed model specifications
#'		\url{http://}
#'
#'@slot kon numeric association rate constant
#'@slot koff numeric dissociation rate constant
#'@slot analyteConcentrations numeric vector for different
#'		analyte concentrations.
#'@slot sd numeric standard deviation of the Gaussian noise. 
#'		by default sd=0, which there is no noise.
#'@slot Rmax numeric the maximum RUs levels at the infinite
#'		high analyte concentration. It is also the maximum level
#'		of ligand immobilized on the chip surface \cr
#'@slot associationLength numeric the time period to run the association
#'@slot dissociationLength numeric the time period to run the dissociation
#'@slot R0 numeric vector the starting RUs for the dissocaiton phase.
#'		optional
#'@slot offset the level shifted at the begining of the dissociation phase
#'		optional
#'@slot Rligand the levels of ligand immobilized on the chip surface. 
#'		optional\cr
#'		If this one is set, then it assumes a variable "Rmax" model. Where 
#'		Rmax for different run is not fixed but rather different due to 
#'		some factors
#'		such as regeration in (Sensiq) or stochasticity (in BioRad Proteon). 
#'		In this model, the immoblized ligand is not fixed from run to run,
#'		but the effienciency or conversion constant is fixed.\cr
#'		This is a vector contains same number of element as the analyteConc.
#'@slot efficiency the conversion constant for a "variable" Rmax model 
#'		optional\cr
#'		see above for explanation.
# #'@seealso 
#'@examples
#'	lgm<-new("LangmuirModel", kon=2E3, koff=0.001, analyteConcentrations=c(1E-6, 2E-6, 1E-5), associationLength=100, 
#'	dissociationLength=100,Rmax=50)
#'	lgm2<-LangmuirModel(1e3, 0.01, analyteConcentrations=c(1e-6, 5e-6, 1e-5), Rmax=30, associationLength=1000, dissociationLength=1000)	
#'	sData<-Simulate(lgm)
#'	#plotting R_AB
#' 	plot(sData[[1]])
#'	#plotting R_A
#'	plot(sData[[2]])
#'	
#'	sData2<-Simulate(lgm2)
#'	#plotting R_AB
#'  plot(sData2[[1]])
#'	#plotting R_A
#'  plot(sData2[[2]])
#' @export
setClass("LangmuirModel",
		representation(kon="numeric",
					   koff="numeric",
					   sd="numeric",
					   analyteConcentrations="numeric", #vector
					   Rmax="numeric", 
					   associationLength="numeric", #option
					   dissociationLength="numeric", #option
					   R0="numeric",#vector optional, if not set the R0 will be getting values from the end of association phase
					   offset="numeric", ###optional, might be used in the future
					   Rligand="numeric", ###optional, only meaningful a variable Rmax model is assumed
					   efficiency="numeric" ###optional, used in the above alternative model
		               ),
		prototype(kon=NULL,
					   koff=NULL,
					   sd=0,
					   analyteConcentrations=numeric(), #vector
					   Rmax=NULL, 
					   associationLength=NULL,
					   dissociationLength=NULL,
					   R0=numeric(),#optional, if not set the R0 will be getting values from the end of association phase
					   offset=0, ###optional, might be used in the future
					   Rligand=numeric(), ###optional, not set by default
					   efficiency=1.0 ###same above, 1.0 by default (100% efficiency)
				)	   
		)
#'@title S4 generic help function to check model validity
#'
#'@description generic function to be called by constructor for 
#'		model validity
#'@details it actually is called by different model method 
#'		setValidty(objects) to check for the validity
#'
#'@param x the model object to be checked for
#'
#' @export 
setGeneric("CheckModelValidity", signature="object",
			function(object) standardGeneric("CheckModelValidity")
			)
#'@describeIn CheckModelValidity check the validity of LangmuirModel
#'
#'@export
setMethod("CheckModelValidity", c("object"="LangmuirModel"),
		function(object)
		{
			#first check for data integrity
				if(is.null(object@kon)){
					cat("kon is not set correctly\n")
					return(FALSE)
				}
				
				if(is.null(object@koff))
				{
					cat("koff is not set correctly\n")
					return(FALSE)
				}
		    #if(is.null(object@sd))
		    #{
		    #  cat("Standard deviation is not set correctly\n")
		    #  return(FALSE)
		    #}
		    if(length(object@analyteConcentrations)<=0)
				{
					cat("analyteConcentrations are empty\n")
					return(FALSE)
				}
				if(sum(object@analyteConcentrations<0)>0)
				{
					cat("analyteConctrations are not set correctly\n")
					return(FALSE)
				}
				if(is.null(object@Rmax))
				{
					cat("Rmax is not set correctly\n")
					return(FALSE)
				}
				if(length(object@R0)>0&&sum(object@R0<0)>0)
				{
					cat("R0 are not set correctly\n")
					return(FALSE)
				}
				if(length(object@R0>0)&&length(object@analyteConcentrations)>0&&length(object@R0)!=length(object@analyteConcentrations))
				{
					cat("R0 and analyteConcentrations are not of equal length\n")
					return(FALSE)
				}
				#check the special case, where R0 has not been set, and we also did run association
				if(is.null(object@associationLength)&&is.null(object@dissociationLength)&&length(object@R0)==0)
				{
					cat("it has specified to run dissociation only, but R0 has not been set")
					return(FALSE)
				}
				#all others are optional.
				
				TRUE
		}
	)
setValidity("LangmuirModel",
			CheckModelValidity
			
		   )
#'@title Constructor for LangmuirModel
#'
#'@seealso \code{\link{LangmuirModel-class}}
#'@examples
#'	lgm<-new("LangmuirModel", kon=2E3, koff=0.001, analyteConcentrations=c(1E-6, 2E-6, 1E-5), associationLength=100, 
#'	dissociationLength=100,Rmax=50)
#'	lgm2<-LangmuirModel(1e3, 0.01, analyteConcentrations=c(1e-6, 5e-6, 1e-5), Rmax=30, associationLength=1000, dissociationLength=1000)	
#'	sData<-Simulate(lgm)
#'	#plotting R_AB
#' 	plot(sData[[1]])
#'	#plotting R_A
#'	plot(sData[[2]])
#'	
#'	sData2<-Simulate(lgm2)
#'	#plotting R_AB
#'  plot(sData2[[1]])
#'	#plotting R_A
#'  plot(sData2[[2]])
#' @export
LangmuirModel<-function(kon,
						koff, sd=0,
						analyteConcentrations,
						Rmax,#="numeric", 
					   associationLength,#="numeric", #option
					   dissociationLength,#="numeric", #option
					   R0=numeric(),#vector optional, if not set the R0 will be getting values from the end of association phase
					   offset=0, ###optional, might be used in the future
					   Rligand=numeric(), ###optional, only meaningful a variable Rmax model is assumed
					   efficiency=1.0 ###optional, used in the above alternative model
					   )
	{
						new("LangmuirModel", kon=kon, koff=koff, sd=sd, analyteConcentrations=analyteConcentrations,
							Rmax=Rmax,
							associationLength=associationLength, dissociationLength=dissociationLength,
							R0=R0, offset=offset,
							Rligand=Rligand, efficiency=efficiency
							)
	}


#doing the simulation
#'@title S4 Geric method to Simluate SPR data
#'
#'@description Geric method to run the simulatoin based on the model
#'
#'@details the each model has different biochemical dynamics to generate
#'		the data. The Langmuir, induced fit and conformational selection
#'		models simulate data analytically, but the two state model simulate
#'		through numerical integration. Please check the detail here
#'		\url{http://}
#'
	
#'@param sampleFreq numeric the time frequency to collect the SPR data
#'@param sd standard deviation used to simulate the nosiy data. 
#'  zero by default for non-noisy data
#'@param fix.ligand currently only used for two different model for 
#'	Rmax: \cr1)model=1, fixed ligand immobilization model. Where Rmax is
#'	kept at a costant level for different run or different channel.
#'	only one Rmax is set through parameter "Rmax". 2)model=2,variable ligand
#'	immobilization model, where Rmax for each run/channel is not fixed
#'	and only "efficincy" is assumed to be constant. \cr if model 
#'	is set to be other values (not 1 or 2), fixed ligand model 
#'	assumed. 
#'@return a list of \code{\link{SensorgramData-class}} data object 
#'		holding the SPR data for different component in the system
#'@seealso 	\code{\link{LangmuirModel-class}}
#'			\code{\link{SensorgramData-class}}
#'@examples
#'	lgm<-new("LangmuirModel", kon=2E3, koff=0.001, analyteConcentrations=c(1E-6, 2E-6, 1E-5), associationLength=100, 
#'	dissociationLength=100,Rmax=50)
#'	lgm2<-LangmuirModel(1e3, 0.01, analyteConcentrations=c(1e-6, 5e-6, 1e-5), Rmax=30, associationLength=1000, dissociationLength=1000)	
#'	sData<-Simulate(lgm)
#'	#plotting R_AB
#' 	plot(sData[[1]])
#'	#plotting R_A
#'	plot(sData[[2]])
#'	
#'	sData2<-Simulate(lgm2)
#'	#plotting R_AB
#'  plot(sData2[[1]])
#'	#plotting R_A
#'  plot(sData2[[2]])
#' @export		   
setGeneric("Simulate", signature="x",
			function(x, sampleFreq=0.01, sd=0, fix.ligand=TRUE,...) standardGeneric("Simulate"))
#'@describeIn Simulate to simulate the SPR data based on a Langmuir Model
setMethod("Simulate", c("x"="LangmuirModel"),#, "sampleFreq"="numeric"),
		function(x,sampleFreq=0.01,sd=0, fix.ligand=TRUE)
		{##running the analytical solution of langmuirModel
			#check the data integrity 
			if(!CheckModelValidity(x))
			{
				return(NULL)
			}
			#check the model and prepare the Rmax
			localRmax<-rep(x@Rmax, length(x@analyteConcentrations))
			if(!fix.ligand)#variable ligand immobilization model
			{
				#check the input for integrity
				if(length(x@Rligand)!=length(x@analyteConcentrations))
				{
					 stop("ERROR: The Rligand has not been set correctly, please check!")
				}
				localRmax<-x@Rligand*x@efficiency
			}
			#else #fixed Rligand model It has been taken care by now
			#{
			#}
			x_sgd<-new("SensorgramData")#for AB
			x_sgd_A<-new("SensorgramData")# for A
			#cat("Rmax:", x@Rmax)
			#x_time<-seq(0,x@associationLength,by=sampleFreq)
			#we also want to check for the validity
			if(x@associationLength>0) #we will do this association phase
			{
			
				x_Ass<-data.frame();
				x_Ass_A<-data.frame();
				x_time<-seq(0,x@associationLength,by=sampleFreq)
				noise<-rnorm(length(x_time),0,sd)
				
				for(i in c(1:length(x@analyteConcentrations)))
				{
					#cat("****i: round",i, "\n")
					
					temp<-localRmax[i]*x@analyteConcentrations[i]/(x@koff/x@kon+x@analyteConcentrations[i])*(1-exp(-1*x_time*(x@kon*x@analyteConcentrations[i]+x@koff)))
					noise<-rnorm(x_time,0,sd)
					#cat("temp:", temp);
					if(i==1)
					{
						x_Ass<-data.frame(Time=x_time, RU1=temp)
						x_Ass_A<-data.frame(Time=x_time, RU1=localRmax[i]-temp)
					}
					else
					{
						x_Ass<-cbind(x_Ass,data.frame(Time=x_time, RU1=temp+noise))
						x_Ass_A<-cbind(x_Ass_A,data.frame(Time=x_time, RU1=localRmax[i]-temp))
					}
				}
				x_sgd@associationData<-x_Ass+noise
				x_sgd_A@associationData<-x_Ass_A
			}#end of association
			#cat("Simulating dissociation phase......\n")
			#print(x_sgd@associationData)
			if(x@dissociationLength>0)#we will do this dissociation phage
			{
				x_Diss<-data.frame();
				x_Diss_A<-data.frame();
				x_time<-seq(0,x@dissociationLength,by=sampleFreq)
				noise<-rnorm(length(x_time),0,sd)
			
								
				for(i in c(1:length(x@analyteConcentrations)))
				{
					x_R0<-0;
					if(length(x@R0)!=0)
					{
						x_R0<-x@R0[i]
					}
					else
					{
						x_R0<-x_sgd@associationData[length(x_sgd@associationData[,1]),i*2]
						#cat("R0:", x_R0,"\n");
					}
					#R02 for RAB_star
					x_R0<-x_R0+x@offset
					
					temp<-x_R0*exp(-1*x_time*x@koff)
					noise<-rnorm(x_time,0,sd)
					
					if(i==1)
					{
#<<<<<<< HEAD
#						x_Diss<-data.frame(Time=x_time, RU1=temp+noise)
#=======
						x_Diss<-data.frame(Time=x_time, RU1=temp)
#>>>>>>> b75c23ddb36900fe1cb6aa401c3aca97e7b7ab0e
						x_Diss_A<-data.frame(Time=x_time, RU1=localRmax[i]-temp)
					}
					else
					{
						x_Diss<-cbind(x_Diss,data.frame(Time=x_time, RU1=temp))
						x_Diss_A<-cbind(x_Diss_A,data.frame(Time=x_time, RU1=localRmax[i]-temp))
					}
				}
				
				x_sgd@dissociationData<-x_Diss+noise
				x_sgd_A@dissociationData<-x_Diss_A
			}
			x_sgd@analyteConcentrations<-x@analyteConcentrations
			x_sgd_A@analyteConcentrations<-x@analyteConcentrations
			return(list(R_AB=x_sgd,R_A=x_sgd_A))
		}#end of function
	)