R/SPRModel_MultiLigand.R

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

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

				
#'@title S4 class of the SPR multi-liganded model
#'
#'@description multi-liganded model object TODO:need to test R0 and fix.ligand=FALSE
#'
#'@details The multi-langed model assumes homogeneous analyte.
#'		It also assumes the ligands do NOT interact with each other.
#'		Each ligand runs independently with the Langmuire model. 
#'		The observed RUs are the sum of all the independent Langmuir
#'		reactions at each time point t.
#'      check the manual for the detailed model specifications
#'		\url{http://}
#'
#'@slot kon numeric vector association rate constants
#'@slot koff numeric vector dissociation rate constants
#'@slot analyteConcentrations numeric vector for different
#'		analyte concentrations.
#'@slot sd numeric vector standard deviation of the Gaussian noise. 
#'		by default sd=0, which there is no noise.
#'@slot Rmax numeric vector 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 numeric vector the level shifted at the begining of the dissociation phase
#'		optional
#'@slot Rligand list of numeric vectors containing the levels of ligand 
#'		immobilized on the chip surface. Each element is one vector of numeric
#'		corresponding to one ligand in the multi-liganded system. Each vector contains
#'		the immobilized levels of this ligand but across different concentrations
#'		of analytes 
#'		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("MultiLigandModel",
		representation(kon="numeric",
					   koff="numeric",
					   sd="numeric",
					   analyteConcentrations="numeric", #vector
					   Rmax="numeric",   associationLength="numeric", #option
					   dissociationLength="numeric", #option
					   R0="list",#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="list", ###optional, only meaningful a variable Rmax model is assumed
					   efficiency="numeric" ###optional, used in the above alternative model
					   #matrix="matrix"
		               ),
		prototype(kon=numeric(),
					   koff=numeric(),
					   sd=0,
					   analyteConcentrations=numeric(), #vector
					   Rmax=numeric(), 
					   associationLength=NULL,
					   dissociationLength=NULL,
					   R0=list(),#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=list(), ###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"="MultiLigandModel"),
		function(object)
		{
			#first check for data integrity
				if(length(object@kon)<=0){
					cat("kon is not set correctly\n")
					return(FALSE)
				}
				if(sum(object@kon<0)>0)
				{
					cat("kons are not set correctly\n")
					return(FALSE)
				}
				if(length(object@koff)<=0)
				{
					cat("koff is not set correctly\n")
					return(FALSE)
				}
				if(sum(object@koff<0)>0)
				{
					cat("koffs are 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(length(object@Rmax)<=0)
				{
					cat("Rmax is not set correctly\n")
					return(FALSE)
				}
			if(sum(object@Rmax<0)>0)
				{
					cat("Rmaxs are not set correctly\n")
					return(FALSE)
				}
			
			#if(missing(object@associationLength)|| missing(object@dissociationLength))
			#	{
			#		cat("association/dissociation running time has not been specified")
			#		return(FALSE)
			#	}
			
			if(object@associationLength<0||object@dissociationLength<0)
			{
				cat("association/dissociation running time has not been specified correct with negative values!")
					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.
			
			#now make sure all the data field are of each length
			#analyteConc is not this range
			lenLigand<-length(object@kon)
			if(lenLigand!=length(object@koff))
			{
				cat("Error:kon and koff are of unequal length!")
				return(FALSE)
			}
			if(length(object@Rmax)>1&&lenLigand!=length(object@Rmax))
			{
				cat("Error:Rmax has not been set to have equal number of element as kon and koff !")
				return(FALSE)
			}
			
			if(class(object@Rligand)!="list")
			{
				cat("Error: Rligand slot has not been specified correctly")
				return(FALSE)
			}
			
			if(class(object@R0)!="list")
			{
				cat("Error: R0 slot has not been specified correctly")
				return(FALSE)
			}
			if(length(object@R0)>0&&length(object@R0)!=length(object@kon))
			{
				cat("Error: Rligand slot has not been specified correctly, inconsist size")
				return(FALSE)
			}
			
			if(is.null(object@associationLength)||object@associationLength<=0)
			{
				stop("Error: associationLength has not been set up correctly!")
			}
			
			if(is.null(object@dissociationLength)||object@dissociationLength<=0)
			{
				stop("Error: dissociationLength has not been set up correctly!")
			}
				TRUE
		}
	)
setValidity("MultiLigandModel",
			CheckModelValidity
			
		   )
#'@title Constructor for LangmuirModel
#'
#'@seealso \code{\link{MultiLigandModel-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
MultiLigandModel<-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=list(), ###optional, only meaningful a variable Rmax model is assumed
					   efficiency=1.0 ###optional, used in the above alternative model
					   )
	{
						new("MultiLigandModel", kon=kon, koff=koff, sd=sd, analyteConcentrations=analyteConcentrations,
							Rmax=Rmax,
							associationLength=associationLength, dissociationLength=dissociationLength,
							R0=R0, offset=offset,
							Rligand=Rligand, efficiency=efficiency
							)
	}


#'@describeIn Simulate to simulate the SPR data based on a MultiLiganded Model
setMethod("Simulate", c("x"="MultiLigandModel"),#, "sampleFreq"="numeric"),
		function(x,sampleFreq=0.01,sd=0, fix.ligand=TRUE)
		{##running the analytical solution of multiligandModel
			#check the data integrity 
			if(!CheckModelValidity(x))
			{
				return(NULL)
			}
			
			#now prepare input Rmax, make a new localRmax,
			#and arrange localRmax as a list,
			#so its elements are vectors, with same number as kon and koff
			#each vectors contains Rmax for each different analyte concentration
			localRmax<-vector("list",length(kon));
			RmaxLigandVec<-x@Rmax
			if(length(RmaxLigandVec)==1)
			{
				# prepare the Rmax
				RmaxLigandVec<-rep(RmaxLigandVec, length(x@kon))
			}
			#expand per each conc
			for(i in 1:length(RmaxLigandVec))
			{
				localRmax[[i]]<-rep(RmaxLigandVec[i],length(x@analyteConcentrations))
			}
			#now we have the a good Rmax list for fix.ligand case,
			
			#try to make a Rmax list when we have a variable ligand case
			if(!fix.ligand)#variable ligand immobilization model
			{
				if(length(x@efficiency)<=0)
				{
					stop("Error: no efficiency has been set")
				}else if(length(x@efficiency)==1)
				{
					x@efficiency<-rep(x@efficiency, length(kon));
				}else if(length(x@efficency)!=length(kon))
				{
					stop("error:efficiency vector is not compatible with other parameters (size)")
				}#else
				#{
				#}
				
				#check the input for integrity
				#remember!!! Rligand is a list with immobilized levels for different concentration as elements
				if(length(x@Rligand)!=length(x@kon))
				{
					 stop("ERROR: The Rligand has not been set correctly, please check!")
				}
				for(i in 1:length(x@Rligand))
				{
					if(is.null(x@Rligand[[i]]))
					{
						stop("Error:the Rligand has not been set correctly, null found!")
					}else if(x@Rligand[[i]]<=1)
					{
						cat("Error: the Rligand has not been set correctly, single value at ",i,"th element")
						return(NULL)
					} else if(length(x@Rligand[[i]])!=x@analyteConcentrations)#&&length(x@Rligand[[i]])>1))
					{
						cat("Error: the Rligand has not been set correctly, incorrect size at ",i,"th element")
						return(NULL)
					} else  #correct case, ignore
					{
						localRmax[[i]]<-x@Rligand[[i]]*x@efficiency[i]	
					}
				}
				#done and we should have set up local Rmax
				#localRmax<-x@Rligand*x@efficiency
			}
			#else #fixed Rligand model It has been taken care by now
			#{
			#}
			
			###checking for R0 input integrity
			localR0<-x@R0
			flag_simulateR0<-TRUE #default to us 
			if(length(x@R0)<=0)
			{
				flag_simulateR0<-TRUE
				localR0<-vector("list", length(x@kon))
				
			}else if(length(x@R0)!=length(x@kon))
			{
				flag_simulateR0<-TRUE
				#discard R0, since R0 is not correctly set
				localR0<-vector("list", length(x@kon))
			}else 
			{
				flag_simulateR0<-FALSE
				#list size is good
				for(i in 1:length(x@R0))
				{
					if(length(x@R0[[i]])!=length(x@analyteConcentrations))
					{
						cat("R0 has not been set correctly, run simulation to get R0")
						flag_simulateR0<-TRUE
						localR0[[i]]<-rep(0,length(x@analyteConcentrations))
					}
				}
			}
			if(flag_simulateR0)
			{
				#list size is good
				for(i in 1:length(localR0))
				{
					
					localR0[[i]]<-rep(0,length(x@analyteConcentrations))
					
				}
			}
			#now R0 is done.
			
			
			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<-0
					tempRmax<-0
					for(j in c(1:length(x@kon)))
					{
						localTemp<-localRmax[[j]][i]*x@analyteConcentrations[i]/(x@koff[j]/x@kon[j]+x@analyteConcentrations[i])*(1-exp(-1*x_time*(x@kon[j]*x@analyteConcentrations[i]+x@koff[j])))
						temp<-temp+localTemp
						tempRmax<-tempRmax+localRmax[[j]][i]
						if(flag_simulateR0)
						{
							localR0[[j]][i]<-localTemp[length(x_time)]
						}
						
					}
					noise<-rnorm(length(x_time),0,sd)
					#cat("temp:", temp);
					if(i==1)
					{
						x_Ass<-data.frame(Time=x_time, RU1=temp+noise)
						x_Ass_A<-data.frame(Time=x_time, RU1=tempRmax-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=tempRmax-temp))
					}
				
				}
				x_sgd@associationData<-x_Ass
				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)))
				{
					temp<-0
					x_R0<-0;
					tempRmax<-0
					for(j in c(1:length(koff)))
					{
					#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<-localR0[[j]][i]+x@offset
					
						temp<-temp+x_R0*exp(-1*x_time*x@koff[[j]])
						noise<-rnorm(length(x_time),0,sd)
						tempRmax<-tempRmax+localRmax[[j]][i]
					}
					if(i==1)
					{
#<<<<<<< HEAD
#						x_Diss<-data.frame(Time=x_time, RU1=temp+noise)
#=======
						x_Diss<-data.frame(Time=x_time, RU1=temp+noise)
#>>>>>>> b75c23ddb36900fe1cb6aa401c3aca97e7b7ab0e
						x_Diss_A<-data.frame(Time=x_time, RU1=tempRmax-temp)
					}
					else
					{
						x_Diss<-cbind(x_Diss,data.frame(Time=x_time, RU1=temp+noise))
						x_Diss_A<-cbind(x_Diss_A,data.frame(Time=x_time, RU1=tempRmax-temp))
					}
				}
				
				x_sgd@dissociationData<-x_Diss
				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
	)