PROGRAM CPF Lifestage
! CPF Lifestage 2015 MR_preg.csl
! Orig from Timchalk 2002
! By C. Timchalk, Battelle PNNL, Richland, WA, USA
! This is based on the final version of the model developed for Dow Agrosciences
! the parameters were experimentally determined or optimized against
! available data sets. See associated data file for needed PROC files
! Updates:
! Metabolic parameters to incorporate scaling
! factors based on experimental age-dependent changes (Smith et al., 2010)
! Amount of tissue esterase (umole) based on
! age-dependent scaling using available published results
! Equation to calc bodywt from age (Hinderliter)
! Incorporated dynamic age-dependent changes to physiology and a generalized 
! framework for human extrapolation (Luecke 2007 & Young 2009)
! Changed QCC flow method to a summation of QC from each compartment.
! GI, diaphragm, dermal and brain metabolism
! Variability estimation added to compartment volumes and enzymology
! Based on Local Sensitivity Analysis results (Poet and Hinderliter)
! Changes to model in January of 2014...
! code was added to allow for dosing of amount via the oral route (vs /kg) (Poet)
! in general, the amount dose variables were written to match the oral route
!also note lines in .m files to setup parameters have been added to zero 
!amount dosed
!Changes in model in January of 2014 to allow for repeat inhalation dosing. (Poet)
!DERMAL and inhalation dosing now controlled by a pulse. Instructions in that area of code (Poet)
!Changes in January of 2015 to include pregnancy (Poet)
!########################################################################

INITIAL

VARIABLE TIME=0.0	! HR
CONSTANT AGE0=0.0	! Initial age (D for rats & Y for humans)
CONSTANT SPECIES=0	! Species dummy variable: 0=male rat, 1=male human
CONSTANT PREG=0.0	!switch to pregnant woman
CONSTANT GDSTART=0.0 !ZERO DAY OF PREGNANCY - DON'T EXCEED 275 (9 MONTHS)

CONSTANT BWSW=0.0 !Switch for age (set =0) or set (set =1) determination of BW
CONSTANT BWST=0.0 !Set value of BW (kg)
! Gompertz equation parameters
CONSTANT B1=0.0;CONSTANT B2=0.0;CONSTANT B3=0.0;CONSTANT B4=0.0;CONSTANT B5=0.0
CONSTANT B6=0.0;CONSTANT B7=0.0;CONSTANT B8=0.0;CONSTANT AGE1=0.0
CONSTANT AGE2=0.0;CONSTANT AGE3=0.0
CONSTANT SEX=0

! Age dependent B4 & B8 for humans (Luecke 2007)
IF ((AGE0 .GE. AGE1) .AND. (AGE0 .LT. AGE2) .AND. (SPECIES .EQ. 1)) THEN
B4x=B4;B8x=0
ELSE IF ((AGE0 .GE. AGE3) .AND. (SPECIES .EQ. 1)) THEN
B4x=0;B8x=B8
ELSE
B4x=0;B8x=0
END IF

! Initial BWT (kg) (Luecke 2007)
IF (AGE0 .GE. AGE2) THEN					 
BWTINIT=B5*EXP((B6/B7)*(1-EXP(-B7*(AGE0-AGE2))))+B8x*(AGE0-AGE3)		
ELSE
BWTINIT=B1*EXP((B2/B3)*(1-EXP(-B3*AGE0)))+B4x*(AGE0-AGE1)
END IF

!INITIALIZE SOME VARIABLES
BWT=(BWTINIT*(1.-BWSW)+VFET+VPLA+VUTR+(VFI-VFINIT)+(VBLI-VBLINIT)+(VSI-VSINIT)+(VRI-VRINIT))+BWST*BWSW
BWTG=(BWTINIT*(1.-BWSW)+BWST*BWSW)*1000			! BWT kg to grams for compartment vol calcs


! Initial Height (HT: fit to CDC growth chart	(human)	
HT=19.994+(5.5408*BWT)-(0.0671*BWT**2)+0.0003*BWT**3


! Initial compartment volumes (frac of BWT)(developed Luecke 2007 & Young 2009; Described in Smith et al., 2013)
VBC=VB0+(VB1*LBWT)+(VB2*LBWT**2)+(VB3*LBWT**3)+(VB4*LBWT**4)+(VB5*LBWT**5)+(VB6*LBWT**6)
								! Brain
VHC=VH0+(VH1*BWTG)+(VH2*BWTG**2)+(VH3*BWTG**3)+(VH4*BWTG**4)	! Liver
VBLC=VBL0+(VBL1*BWTG)+(VBL2*BWTG**2)				! Blood
VFC=VF0+(VF1*BWTG)+(VF2*BWTG**2)+(VF3*BWTG**3)+(VF4*BWTG**4)+(VF5*BWTG**5)+(VF6*BWTG**6)	
VAC=VA0+(VA1*BWTG)+(VA2*BWTG**2)+(VA3*BWTG**3)+(VA4*BWTG**4)+(VA5*BWTG**5)+(VA6*BWTG**6)	
								! Adipose
VKC=VK0+(VK1*BWTG)+(VK2*BWTG**2)				! Kidney (in rapid)
VSPC=VSP0+(VSP1*BWTG)+(VSP2*BWTG**2) 				! Spleen (in rapid)
VLC=VL0+(VL1*BWTG)+(VL2*BWTG**2)				! Lung (in rapid)
VGIC=VGI0!+(VGI1*BWTG)+(VGI2*BWTG**2)+(VGI3*BWTG**3)		! GI, both stomach and intestines (in rapid)
VMC=VM0+(VM1*BWTG)+(VM2*BWTG**2)+(VM3*BWTG**3)+(VM4*BWTG**4)			! Muscle (in slow)
VSKC=VSK0+(VSK1*BWTG)+(VSK2*BWTG**2)+(VSK3*BWTG**3)+(VSK4*BWTG**4)+(VSK5*BWTG**5)
								! Skin (in slow) note - for dermal this is "external" compartment
CONSTANT FATVAR = 0.0
CONSTANT VOLUNCERT = 0.0
CONSTANT KUNCERT= 0.0
CONSTANT VHCv=0.0
CONSTANT VBLv=0.0
QHCV=0.0
VFCv=GAUSS(0.,FATVAR*VFC)
VFCv = RSW(VFC+VFCv .LT. 0.1*VFC, 0.1*VFC, VFCv) ! check to ensure that vol is > 0
CONSTANT VBMC=0.0						! Bone marrow (in slow)
VRC=VKC+VSPC+VLC+VGIC				! Rapid
VSC=VMC+VSKC+VBMC+(VAC-VFC)					! Slow
CONSTANT VDC=0.0						! Diaphram
CONSTANT VPC=0.0						! Pancreas (in rapid)
CONSTANT VINTC=0.0						! Small intestine, for metabolism only

! Volume variability, set VOLUNCERT=0 to turn off
VSCu=GAUSS(0.,VOLUNCERT*VSC)
VRCu=GAUSS(0.,VOLUNCERT*VRC)
VBCu=GAUSS(0.,VOLUNCERT*VSC)
KPLU=GAUSS(0.,KUNCERT*KPL)
!VHCu=GAUSS(0.,VOLUNCERT*VHC) !This can be used if total uncertainty is desired, or the next line to vary just hepatic volume based on P3M
VHCu=GAUSS(0.,VHCv*VHC)
VFCu=GAUSS(0.,VOLUNCERT*VFC)
VBLCu=GAUSS(0.,VBLV*VBLC)
VDCu=GAUSS(0.,VOLUNCERT*VDC)
VINTCu=GAUSS(0.,VOLUNCERT*VINTC)
VSCu = RSW(VSC+VSCu .LT. 0.1*VSC, 0.1*VSC, VSCu) ! check to ensure that vol is > 0
VRCu = RSW(VRC+VRCu .LT. 0.1*VRC, 0.1*VRC, VRCu) ! check to ensure that vol is > 0
VBCu = RSW(VBC+VBCu .LT. 0.1*VBC, 0.1*VBC, VBCu) ! check to ensure that vol is > 0
VHCu = RSW(VHC+VHCu .LT. 0.1*VHC, 0.1*VHC, VHCu) ! check to ensure that vol is > 0
VFCu = RSW(VFC+VFCu .LT. 0.1*VFC, 0.1*VFC, VFCu) ! check to ensure that vol is > 0
VBLCu = RSW(VBLC+VBLCu .LT. 0.1*VBLC, 0.1*VBLC, VBLCu) ! check to ensure that vol is > 0
VDCu = RSW(VDC+VDCu .LT. 0.1*VDC, 0.1*VDC, VDCu) ! check to ensure that vol is > 0
VINTCu = RSW(VINTC+VINTCu .LT. 0.1*VINTC, 0.1*VINTC, VINTCu) ! check to ensure that vol is > 0
KPLU = RSW(KPL+KPLu .LT. 0.1*KPL, 0.1*KPL, KPLu) ! check to ensure > 0

VSC=VSC+VSCu
VRC=VRC+VRCu
VBC=VBC+VBCu
VHC=VHC+VHCu
VFC=VFC+VFCu
VBLC=VBLC+VBLCu
VDC=VDC+VDCu
VINTC=VINTC+VINTCu
KPL=KPL+KPLU


! Compartment volume polynomial parameters
CONSTANT VB0=0.0;CONSTANT VB1=0.0;CONSTANT VB2=0.0;CONSTANT VB3=0.0;CONSTANT VB4=0.0
CONSTANT VB5=0.0;CONSTANT VB6=0.0
CONSTANT VH0=0.0;CONSTANT VH1=0.0;CONSTANT VH2=0.0;CONSTANT VH3=0.0;CONSTANT VH4=0.0
CONSTANT VBL0=0.0;CONSTANT VBL1=0.0;CONSTANT VBL2=0.0
CONSTANT VF0=0.0;CONSTANT VF1=0.0;CONSTANT VF2=0.0;CONSTANT VF3=0.0;CONSTANT VF4=0.0
CONSTANT VF5=0.0;CONSTANT VF6=0.0
CONSTANT VA0=0.0;CONSTANT VA1=0.0;CONSTANT VA2=0.0;CONSTANT VA3=0.0;CONSTANT VA4=0.0
CONSTANT VA5=0.0;CONSTANT VA6=0.0
CONSTANT VK0=0.0;CONSTANT VK1=0.0; CONSTANT VK2=0.0
CONSTANT VSP0=0.0;CONSTANT VSP1=0.0;CONSTANT VSP2=0.0
CONSTANT VL0=0.0;CONSTANT VL1=0.0;CONSTANT VL2=0.0
CONSTANT VGI0=0.0;CONSTANT VGI1=0.0;CONSTANT VGI2=0.0;CONSTANT VGI3=0.0
CONSTANT VM0=0.0;CONSTANT VM1=0.0;CONSTANT VM2=0.0;CONSTANT VM3=0.0;CONSTANT VM4=0.0
CONSTANT VSK0=0.0;CONSTANT VSK1=0.0;CONSTANT VSK2=0.0;CONSTANT VSK3=0.0;
CONSTANT VSK4=0.0;CONSTANT VSK5=0.0

! Compartment volumes (L)
VB=VBC*BWT
VH=VHC*BWT
VBLN=VBLC*BWT		!For pregnancy, blood volume is approximated by the combination of plasma and RBC, which increase over GD
								!VBLN is used to set to age-specific blood volume to start
VAB=0.46*VBL
VRN=VRC*BWTINIT
VSN=VSC*BWTINIT
VD=VDC*BWT
VINT=VINTC*BWT
VLU=VLC*BWT
VFN=VFC*BWT !For pregnancy, fat increases over GD, VFN is used to set to the age-specific fat volume to start

IF (SPECIES .EQ. 0) THEN	! Log of BWT for rat brain calc
LBWT=LOG(BWT)			
ELSE
LBWT=BWT*1000
END IF

! Blood flow to each compartment (L/h/kg) 
CONSTANT QFC=0.0	! To fat
CONSTANT QHC=0.0	! To liver
CONSTANT QSC=0.0 	! To slow
CONSTANT QBC=0.0	! To brain
CONSTANT QDC=0.0   	! To diaphram
CONSTANT QRC=0.0	! To rapid

!Variability to to activity level
CONSTANT WORK= 0.0 !Work in Watts
DVO2=WORK*10.1/1000.;
! Blood Flow adjustments
QFCv	=	0.017*DVO2*60.
!QHCv	=	-0.19*DVO2*60.  !FOR WORK, OR FOR VARIABILITY COMMENT OUT and use following line
QHCvv=GAUSS(0.,QHCv*QHC)
QSCv  =  7.38*DVO2*60.
QRCv  =  -0.08*DVO2*60.

! VOLUME OF EXPOSED SKIN - MODEL 1
VSKCC=VSKC*(SA/TSA)	! FRACTION SKIN TISSUE EXPOSED 

! PART COEF CHLORPYRIFOS (C) MODEL 1 (POULIN, KRISHNAN; 1995)
CONSTANT PFCN=0.0	! FAT/BLOOD
CONSTANT PHC=0.0	! LIVER/BLOOD
CONSTANT PRC=0.0	! RICHLY PERFUSED/BLOOD
CONSTANT PSC=0.0	! SLOWLY PERFUSED/BLOOD
CONSTANT PDC=0.0	! DIAPHRAM/BLOOD
CONSTANT PBC=0.0	! BRAIN/BLOOD
CONSTANT PBL=0.0	! BLOOD/BLOOD
CONSTANT PLU=0.0	! LUNG/BLOOD
CONSTANT PLAIR=2000	!NOT VOLATILE
CONSTANT PSKL=0.0  ! THIS IS FOR EXPOSURE ONLY AND REPRESENTS QSKC/PSKC
CONSTANT PPLA=0.0!Placenta
CONSTANT PAFC=1		 !FETAL DIFFUSION - THIS IS 0.01 IN GENTRY MODEL, 19 IN lOWE ET AL 2010

! PROTEIN BINDING CHLORPYRIFOS (C) MODEL 1 (Lowe 2009)
CONSTANT FBC=0.0	! FRAC BOUND CPF IN BLOOD
CONSTANT FBO=0.0	! FRAC BOUND CPF-oxon in blood 
CONSTANT FUBO=0.0  !FRAC UNBOUND
CONSTANT FUBC=0.0
CONSTANT HCTN=0.0	! HEMATOCRIT FRACTION RBC

! PART COEF OXON (O) MODEL 2 (Lowe 2009)
CONSTANT PFON=0.0	! FAT/BLOOD
CONSTANT PHO=0.0	! LIVER/BLOOD
CONSTANT PRO=0.0	! RICHLY PERFUSED/BLOOD
CONSTANT PSO=0.0	! SLOWLY PERFUSED/BLOOD
CONSTANT PDO=0.0	! DIAPHRAM/BLOOD
CONSTANT PBO=0.0	! BRAIN/BLOOD
CONSTANT PBLO=0.0	! BLOOD/BLOOD
CONSTANT PLUO=0.0
CONSTANT PLAIRO=2000
CONSTANT PSKO=0.0
CONSTANT PPLAO=0.0



! Metabolic parameters (Smith 2011, Poet 2004)
CONSTANT KMHHCO=0.0	! Hepatic Km CPF -> CPF-oxon (umol/L)
CONSTANT VMHHCO=0.0	! Hepatic Vmaxc CPF -> CPF-oxon (umol/hr/kg tissue)
CONSTANT KMHCP=0.0	! Hepatic Km CPF -> TCPy (umol/L)
CONSTANT VMHCP=0.0	! Hepatic Vmaxc CPF -> TCPy (umol/hr/kg tissue)
CONSTANT KMLST=0.0	! Hepatic Km CPF-oxon -> TCPy (umol/L)
CONSTANT VML=0.0	! Hepatic Vmaxc CPF-oxon -> TCPy (umol/hr/kg tissue)

CONSTANT VMBL1=0.0;CONSTANT VMBL2=0.0;CONSTANT VMBL3=0.0 
VMBL=VMBL1/(1+EXP((VMBL2-AGE)/VMBL3))		! Blood Vmaxc CPF-oxon -> TCPy (umol/hr/L plasma), Dynamic: logistic func MOVED FOR PREG
!CONSTANT VMBL=430000  !Comment out line above AND 772 to run bootstrap wtih assigned VMBL values)
CONSTANT KMBLST=0.0	! Blood Km CPF-oxon -> TCPy (umol/L)

CONSTANT VINTOC=0.0	! Small intestine Vmaxc of CPF -> CPF-oxon (umol/hr/kg liver)
CONSTANT KMINTO=0.0	! Small intestine Km of CPF -> CPF-oxon (umol/L)
CONSTANT VINTTC=0.0	! Small intestine Vmaxc of CPF -> TCPy (umol/hr/kg tissue)
CONSTANT KMINTC=0.0	! Small intestine Km of CPF -> TCPy (umol/L)
CONSTANT VINTOXC=0.0	! Small intestine Vmaxc of CPF-oxon -> TCPy (umol/hr/kg tissue)
CONSTANT KMINTOX=0.0	! Small intestine Km of CPF-oxon -> TCPy (umol/L)

CONSTANT VMCBCO=0.0	! Brain Vmaxc of CPF -> CPF-oxon (umol/h/kg)
CONSTANT KMBCO=0.0	! Brain Km of CPF -> CPF-oxon (umol/L)
CONSTANT VMCBCP=0.0	! Brain Vmaxc of CPF -> TCPy (umol/h/kg)
CONSTANT KMBCP=0.0	! Brain Km of CPF -> TCPy (umol/L)
CONSTANT KMBCPT=0.0	!Brain oxon ->TCPy
CONSTANT VMBCPTC=0.0 !Brain oxon ->TCPy

CONSTANT VMAXOLUC=0.0	!LUNG VMAXC CPF->OXON
CONSTANT KMOLU=0.0	!LUNG CPF->OXON
CONSTANT VMAXTLUC=0.0	!LUNG VMAXC CPF->TCP
CONSTANT KMTLU=0.0	!LUNG CPF->TCP
CONSTANT VMAXPLUC=0.0	!LUNG VMAXC OXON->TCP
CONSTANT KMPLUC=0.0	!LUNG OXON->TCP
CONSTANT KMUC=0.87    !kMC hr-1 FURUBAYASHI et al Biol. Pharm. Bull. 2007 = 0.0145 min-1

!diaphragm metabolism parameters, Vmaxc (umol/hr/g tissue), Km (uM)copied from brain
CONSTANT VMCDPCO=0.0		! CPF -> CPF-oxon
CONSTANT KMDPCO=0.0		! Fitted to brain ChEA, 1/10 hepatic
CONSTANT VMCDPCP=0.0		! CPF -> TCPy, Extrap from Layshock 2009 unpub
CONSTANT KMDPCP=0.0		! Same as hepatic

!dermal metabolism parameters, Vmaxc (umol/hr/g tissue), Km (uM)copied from brain
CONSTANT VMCDERMCO=0.0		! CPF -> CPF-oxon
CONSTANT KMDERMCO=0.0		! Fitted to brain ChEA, 1/10 hepatic
CONSTANT VMCDERMCP=0.0		! CPF -> TCPy, Extrap from Layshock 2009 unpub
CONSTANT KMDERMCP=0.0		! Same as hepatic
CONSTANT KMDERMEST=0.0	! Km CPF-oxon -> TCPy (umol/L), same as hepatic
CONSTANT VMDERMESTC=0.0	! dermal Vmaxc CPF-oxon -> TCPy (umol/hr/kg tissue) 

! Scaling Vmaxc to Vmax (umol/h) 
VMHHON=VMHHCO*VH		     	! Hepatic CPF -> CPF-oxon
VMHCN=VMHCP*VH			       ! Hepatic CPF -> TCPy
VMLSTN=VML*VH
VMBLSTN=VMBL*VBL*(1-HCT)
!VMLST=VMLN*VH		        	! Hepatic CPF-oxon -> TCPy moved to dynamic for pregnancy control
VINTO=VINTOC*VINT		   	! Sm int. CPF -> CPF-oxon
VINTT=VINTTC*VINT			   ! Sm int. CPF -> TCPy
VINTOX=VINTOXC*VINT			 ! Sm int. CPF-oxon -> TCPy
VMBCO=VMCBCO*VB			     ! Brain CPF -> CPF-oxon
VMBCP=VMCBCP*VB			     ! Brain CPF -> TCPy
VMBCPT=VMBCPTC*VB
VMAXOLU=VMAXOLUC*VLU		! LUNG CPF->OXON
VMAXTLU=VMAXTLUC*VLU		! LUNG CPF->TCP
VMAXPLU=VMAXPLUC*VLU		!LUNG OXON->TCP
VMCDPO=VMCDPCO*VD	!Diaphragm CPF-> oxon
VMCDPP=VMCDPCP*VD !Diaphragm CPF-> TCPy
VMDERMO=VMCDERMCO*VSKCC		     	! DERMAL CPF -> CPF-oxon
VMDERMC=VMCDERMCP*VSKCC						! DERMAL CPF -> TCPy
VMDERMEST=VMDERMESTC*VSKCC			! DERMAL ESTERASE OXON -> TCPy


! Enzyme turnover rates (umol hydrolysis/h/umol active sites)
CONSTANT TRCE=0.0		! ACHE
CONSTANT TRBE=0.0		! BCHE
CONSTANT TRCR=0.0		! CAE

! Enzyme activities (umol/hr/kg tissue)
CONSTANT BACHE=0.0		! BRAIN AChE
CONSTANT DACHE=0.0		! DIAPHRAGM AChE
CONSTANT HACHE=0.0		! HEPATIC AChE
CONSTANT BLACHE=0.0		! Plasma AChE
CONSTANT RBCHE=0.0		! RBC AChE
CONSTANT LUACHE=0.0		! LUNG AChE

CONSTANT BRCE=0.0    ! Brain CaE
CONSTANT DACE=0.0		 ! DIAPHRAGM CaE
CONSTANT HECE=0.0    ! Hepatic CaE
CONSTANT PLOCE=0.0		! Plasma CaE
CONSTANT LUECE=0.0    ! LUNG CaE

CONSTANT BBUCE=0.0		! BRAIN BuChE
CONSTANT DBUCE=0.0		! DIAPHRAGM BuChE
CONSTANT HBUCE=0.0		! LIVER BuChE
CONSTANT BLBUCE=0.0		! Plasma BuChE
CONSTANT LUBUCE=0.0		!LUNG BuChE

! HEPATIC B-CHOLINE (C) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDHCE=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIHCE=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRHCE=0.0	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KAHCE=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! HEPATIC B-BUTYRL (B) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDHBE=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIHBE=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRHBE=0.0	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KAHBE=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! HEPATIC B-CARBOXYL (CR) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDHCR=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIHCR=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRHCR=0.0	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KAHCR=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! BRAIN B-CHOLINE (C) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDBCE=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIBCE=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRBCE=0.0	! HR-1; REACTIVATION       OF ESTERASE
CONSTANT KABCE=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! BRAIN B-BUTYRL (B) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDBBE=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIBBE=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRBBE=0.0	! HR-1;REACTIVATION OF ESTERASE
CONSTANT KABBE=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! BRAIN B-CARBOXYL (CR) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDBCR=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIBCR=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRBCR=0.0	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KABCR=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! LUNG B-CHOLINE (C) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDLUCE=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KILUCE=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRLUCE=0.0	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KALUCE=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! LUNG B-BUTYRL (B) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDLUBE=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KILUBE=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRLUBE=0.0	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KALUBE=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! LUNG B-CARBOXYL (CR) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDLUCR=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KILUCR=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRLUCR=0.0	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KALUCR=0.0	! HR-1; RATE FOR AGING OF ESTERASE


! DIAPHRAM B-CHOLINE (C) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDDCE=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIDCE=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRDCE=0.0	! HR-1;REACTIVATION OF ESTERASE
CONSTANT KADCE=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! DIAPHRAM B-BUTYRL (B) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDDBE=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIDBE=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRDBE=0.0	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KADBE=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! DIAPHRAM B-CARBOXYL (CR) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDDCR=0.0	! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIDCR=0.0	! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRDCR=0.0	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KADCR=0.0	! HR-1; RATE FOR AGING OF ESTERASE

! BLOOD/PLASMA B-CHOLINE (C) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDBLCE=0.0		! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIBLCE=0.0		! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRBLCE=0.0		! HR-1; REACTIVATION OF ESTERASE
CONSTANT KABLCE=0.0		! HR-1; RATE FOR AGING OF ESTERASE

! BLOOD/PLASMA B-BUTYRL (B) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDBLBE=0.0		! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIBLBE=0.0		! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRBLBE=0.0	 	! HR-1; REACTIVATION OF ESTERASE
CONSTANT KABLBE=0.0 		! HR-1; RATE FOR AGING OF ESTERASE

! BLOOD/PLASMA B-CARBOXYL (CR) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDBLCR=0.0		! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIBLCR=0.0		! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRBLCR=0.0		! HR-1; REACTIVATION OF ESTERASE
CONSTANT KABLCR=0.0		! HR-1; RATE FOR AGING OF ESTERASE

! RBC B-CHOLINE (C) ESTERASE RATE CONSTANTS - MODEL 2
CONSTANT KDRBCE=0.0		! HR-1; DEGRADATION OF ESTERASE
CONSTANT KIRBCE=0.0		! UMOLES-1 L HR-1; INHIBITION OF ESTERASE
CONSTANT KRRBCE=0.0		! HR-1; REACTIVATION OF ESTERASE
CONSTANT KARBCE=0.0		! HR-1; RATE FOR AGING OF ESTERASE

! Enzyme activity (umol/h)
! AChE
SBACH=BACHE*VB
SDACH=DACHE*VD
SHACH=HACHE*VH	
SBLACH=BLACHE*VBL*(1-HCT)
SRBACH=RBCHE*VBL*HCT
SLUACH=LUACHE*VLU
! CAE
SBRCE=BRCE*VB	
SDACE=DACE*VD		
SHECE=HECE*VH			
SPLOCE=PLOCE*VBL*(1-HCT)
SLUECE=LUECE*VLU

! BuChE
SBBUC=BBUCE*VB
SDBUC=DBUCE*VD	
SHBUC=HBUCE*VH	
SBLBUC=BLBUCE*VBL*(1-HCT)
SLUBUC=LUBUCE*VLU

! CALCULATION OF ESTERASE BINDING SITES (UMOLES) (MAXWELL, 1987)
IBCE=SBACH/TRCE		! TOTAL BRAIN ACHE (UMOLES)
IDCE=SDACH/TRCE		! TOTAL DIAPHRAGM ACHE (UMOLES)
IHCE=SHACH/TRCE		! TOTAL HEPATIC ACHE (UMOLES)
IBLCE=SBLACH/TRCE	! TOTAL PLASMA ACHE (UMOLES)
IRBCE=SRBACH/TRCE	! Total RBC AChE (umol)
ILUCE=SLUACH/TRCE		! TOTAL LUNG ACHE (UMOLES)

IBBE=SBBUC/TRBE		! TOTAL BRAIN BUCHE (UMOLES)
IDBE=SDBUC/TRBE		! TOTAL DIAPHRAGM BUCHE (UMOLES)
IHBE=SHBUC/TRBE		! TOTAL HEPATIC BUCHE (UMOLES)
IBLBE=SBLBUC/TRBE	! TOTAL PLASMA BUCHE (UMOLES)
ILUBE=SLUBUC/TRBE		! TOTAL LUNG BUCHE (UMOLES)

IBCR=SBRCE/TRCR		! TOTAL BRAIN CAE (UMOLES)
IDCR=SDACE/TRCR		! TOTAL DIAPHRAGM CAE (UMOLES)
IHCR=SHECE/TRCR		! TOTAL HEPATIC CAE (UMOLES)
IBLCR=SPLOCE/TRCR	! TOTAL PLASMA CAE (UMOLES)
ILUCR=SLUECE/TRCR		! TOTAL LUNG CAE (UMOLES)


! CALCULATION OF ENZ SYN RATE (UMOLE/HR)
! LIVER
KSHCE=IHCE*KDHCE	! SYN RATE OF HEPATIC ACHE (UMOLE/HR)
KSHBE=IHBE*KDHBE	! SYN RATE OF HEPATIC BCHE (UMOLE/HR)
KSHCR=IHCR*KDHCR	! SYN RATE OF HEPATIC BCHE (UMOLE/HR)
! BRAIN
KSBCE=IBCE*KDBCE	! SYN RATE OF BRAIN ACHE (UMOLE/HR)
KSBBE=IBBE*KDBBE	! SYN RATE OF BRAIN BCHE (UMOLE/HR)
KSBCR=IBCR*KDBCR	! SYN RATE OF BRAIN BCHE (UMOLE/HR)
! LUNG
KSLUCE=ILUCE*KDLUCE	! SYN RATE OF LUNG ACHE (UMOLE/HR)
KSLUBE=ILUBE*KDLUBE	! SYN RATE OF LUNG BCHE (UMOLE/HR)
KSLUCR=ILUCR*KDLUCR	! SYN RATE OF LUNG BCHE (UMOLE/HR)
! DIAPHRAGM
KSDCE=IDCE*KDDCE	! SYN RATE OF DIAPHRAGM ACHE (UMOLE/HR)
KSDBE=IDBE*KDDBE	! SYN RATE OF DIAPHRAGM BCHE (UMOLE/HR)
KSDCR=IDCR*KDDCR	! SYN RATE OF DIAPHRAGM BCHE (UMOLE/HR)
! PLASMA
KSBLCE=IBLCE*KDBLCE	! SYN RATE OF PLASMA ACHE (UMOLE/HR)
KSBLBE=IBLBE*KDBLBE	! SYN RATE OF PLASMA BCHE (UMOLE/HR)
KSBLCR=IBLCR*KDBLCR	! SYN RATE OF PLASMA BCHE (UMOLE/HR)
! RBC
KSRBCE=KDRBCE*IRBCE	! SYNTHESIS RATE OF RBC (umol/h)


! ACUTE BOLUS ORAL DOSING OF CPF - MODEL 1
CONSTANT ORALC=0.0	! CPF single oral bolus dose (ug/kg)
CONSTANT ORALP=0.0	! TCPy oral bolus dose (ug/kg)
CONSTANT ORALO=0.0	! CPF-oxon oral bolus dose (ug/kg)
CONSTANT KAI=0.0	! 1ST ORDER ABSORPTION RATE INTESTINE HR-1
CONSTANT KSI=0.0	! 1ST ORDER TRANSFER STOM TO INTEST HR-1
CONSTANT KAS=0.0   !RAT ONLY, STOM TO LIVER HR-1
CONSTANT KIHP=0.0	! Transfer of TCPy from sm. intest to compartment 1 (/h)
CONSTANT KOS=0.0	! Transfer of CPF-oxon from sm. intest to liver (/h)
CONSTANT KOSI=0.0 !Transfer of CPF-oxon (from an oral dose) to the intestines, note in CPF model for ease
CONSTANT FA=0.0		! FRACTIONAL ABSORPTION FROM GI TRACT


! MOLECULAR WEIGHTS
CONSTANT MWC=350.6	! MW FOR CHLORPYRIFOS - MODEL 1 (G/MOL)
CONSTANT MWO=334.5	! MW FOR CHLORPYRIFOS OXON  - MODEL 1 & 2 (G/MOL)
CONSTANT MWP=198.4	! MW FOR TCP - MODEL 1 & 2 (G/MOL)

! Oral Dosing - dose an undefined number of times
! Put times in the array DT (in hours) (absolute time, NOT time since last dose)
! Put the oral dose of chlorpyrifos in the array DORAL (UG/KG)(DORALO FOR OXON)
!for dosing AMOUNT instead of /kg (ug) put the oral dose of CPF (AMTCDOSE) OR OXON (AMTODOSE)
! Add one additional dose of zero at a time greater than TSTOP to keep the model running
DIMENSION DT(2000), DORAL(2000), DSC(2000), DORALO(2000), AMTCDOSE(2000), AMTODOSE(2000)
INTEGER indx
indx=1  			! indx is a counter used to track multiple doses
SCHEDULE DOSEX .AT. DT(indx)
ODOSE=0.0
SCDOSE=0.0
ODOSEO=0.0
ADOSEC=0.0
ADOSEO=0.0



! MOLAR CONVERSION & FRACTIONAL ABSORPTION MODIFICATION
ODOSEC=ORALC*(BWT/MWC)*FA	! CPF dose conversion (umol)
ODOSEP=ORALP*(BWT/MWP)*FA	! TCPy dose coversion (umol)
ODOSEO=ORALO*(BWT/MWO)*FA	! CPF-oxon dose conversion (umol)

! Oral infusion AKA FEED
CONSTANT FEED=0.0		! Feed dose (ug/kg)
CONSTANT FEEDR=0.0 		! Rate of feed intake, est @ 1.3 day-1 (Saghir 2006)
CONSTANT FEEDT=0.0 		! TIME OF FEEDING
FEEDOSE=FEED*FEEDR*BWT/24 		! ug/day
FEDOSE1=FEEDOSE/MWC		! UMOLES/day - NO FA, ASSUME 100%
COMPDOSE=FEEDOSE/BWT*FEEDT 	! TOTAL UG/KG DOSED
FEEDBALANCE=COMPDOSE*(BWT/MWC)	! CPF dose conversion (umol)

! IV DOSE INFUSION
CONSTANT IV=0.0			! IV DOSE (uG/KG BWT)			
CONSTANT IVS=0.0		! IV START TIME (HR)
CONSTANT IVE=0.01		! IV END TIME (HR)
IVD=IVE-IVS			! IV INFUSION DURATION (HR)			
TIV=IV*BWT/MWC			! TOTAL CPF INJECTED (uMOL)
IVR=TIV/IVD			! RATE OF CPF DOSE INFUSION (uMOL/HR)


! SUBCUTANEOUS DOSE (JNS)
CONSTANT SCDC=0.0		! SUBQ DOSE (uG/KG BWT)
CONSTANT KSC=0.0		! 1ST ORDER TRANSFER RATE (HR^-1)
CONSTANT KSP=0.0		! 1ST ORDER TRANSFER RATE FROM SUBQ TO PERIPHERAL SUBQ
CONSTANT KPS=0.0		! 1ST ORDER TRANSFER RATE FROM PERIPHERAL SUBQ TO SUBQ
DSCDC=((SCDC*BWT)/MWC)		! AMOUNT OF CPF DOSE (uMOL)

!DERMAL DOSE (TSP)
CONSTANT DERM=0.0	!UG/KG or ug, depending on which line is used below 549 or 550
CONSTANT SAF=0.0   !FRACTION TOTAL SKIN EXPOSED
CONSTANT VLIQ=0.0	!L
CONSTANT KPL=0.0	!cm/hr
CONSTANT KDIS = 0. ! 1st order dissipation loss from skin to environment, 1/day
CONSTANT FAD=0.0	!FRACT DERM ABSORBED, should be 100% available = 1

!REPIRATORY DOSE
CONSTANT CONCPPM = 0.0    ! Inhaled conc (ppm)
CONSTANT CONCMGM = 0.0		!Inhaled conc (mg/m3)
CONSTANT CONCPPMO=0.0
CONSTANT CONCMGMO=0.0
CONSTANT CCON = 0.0    ! For closed chamber
CONSTANT CCOFF = 1.0    ! For closed chamber
CONSTANT VCHCC = 6400.0 ! Closed chamber volume (mL)
VCHA= VCHCC-N*BWT	! ...corrected for animal volume
CONSTANT N = 0      ! Number of animals per chamber
CONSTANT KLOSS = 0.0 ! Loss from chamber, 1/min
CONSTANT DAYSWK = 1.0    ! Number of exposure days per week
CONSTANT TMAX = 1440.0 ! Maximum time for exposures
CONSTANT FRACIN=0.0
CONSTANT s2=0											!INHALATION ON
CONSTANT p2=3											!DERMAL EXPOSURE

CONSTANT D3=5  !5 DAYS/WEEK FOR DERMAL EXPOSURE
CONSTANT YEAR=364  !DAYS! TO STOP DERMAL and INHALTION EXPOSURE AT SOME POINT DURING A YEAR (THIS COULD ALSO BE MONTHS...)
CONSTANT W2=21  ! DAYS! THIS COMBINES WITH YEAR - 21 HR=3 WEEKS OF EXPOSURE.
CONSTANT D3IN=5  !5 DAYS/WEEK FOR inhalation EXPOSURE
CONSTANT W2IN=21  ! THIS COMBINES WITH YEAR - 21 WEEKS OF EXPOSURE.- INHALATION
CONSTANT P2IN=3	!INHALATION HR/DAY
CONSTANT S2IN=0 !TO OFFSET DSOING START, PROBABLY NEVER SET TO ANYTHING BUT 0
CONSTANT VVOL=0.06  !URNARY OUTPUT, L/HR -Hayes et al, 2015 from NHANES data



! COMPARTMENTAL MODEL FOR TCP ELIMINATION
CONSTANT KE=0.0		! ELIMINATION RATE (HR-1)
V1=0.176*BWT^1.01	! TCPy vol of distribution (L)  0.176*bwt^1.01 - see Smith et al, 2013


! Enzymatic variability
CONSTANT VMLv=0.0
CONSTANT VMBLv=0.0
CONSTANT VMHCPv=0.0
CONSTANT VMHHCOv=0.0
CONSTANT VMCBCPv=0.0
CONSTANT VMCBCOv=0.0
CONSTANT VINTTv=0.0
CONSTANT VINTOv=0.0

! lognormal distributions
VMLvv=exp(GAUSS(LOG(VMLSTN)-(SQRT(LOG((VMLv)**2.+1))**2.)/2.,SQRT(LOG((VMLv)**2.+1))))-VMLSTN
VMBLvv=((VMLvv+VMLST)/VMLSTN*VMBLSTN+GAUSS(0.,VMBLv))-VMBLSTN
VMHCPvv=exp(GAUSS(LOG(VMHCN)-(SQRT(LOG((VMHCPv)**2.+1))**2.)/2.,SQRT(LOG((VMHCPv)**2.+1))))-VMHCN
VMHHCOvv=exp(GAUSS(LOG(VMHHON)-(SQRT(LOG((VMHHCOv)**2.+1))**2.)/2.,SQRT(LOG((VMHHCOv)**2.+1))))-VMHHON
VMCBCPvv=exp(GAUSS(LOG(VMBCP)-(SQRT(LOG((VMCBCPv)**2.+1))**2.)/2.,SQRT(LOG((VMCBCPv)**2.+1))))-VMBCP
VMCBCOvv=exp(GAUSS(LOG(VMBCO)-(SQRT(LOG((VMCBCOv)**2.+1))**2.)/2.,SQRT(LOG((VMCBCOv)**2.+1))))-VMBCO
VINTOvv=exp(GAUSS(LOG(VINTO)-(SQRT(LOG((VINTOv)**2.+1))**2.)/2.,SQRT(LOG((VINTOv)**2.+1))))-VINTO
VINTTvv=exp(GAUSS(LOG(VINTT)-(SQRT(LOG((VINTTv)**2.+1))**2.)/2.,SQRT(LOG((VINTTv)**2.+1))))-VINTT
   
!If the Vmax is less than 1% of the nominal, set equal to 1% of nonimal (by setting error to 99% of nominal)
VMLvv = RSW(VMLvv+VMLSTN .LT. .01*VMLSTN, -0.99*VMLSTN, VMLvv)
VMBLvv = RSW(VMBLvv+VMBLSTN .LT. .1*VMBLSTN, -0.90*VMBLSTN, VMBLvv)		!NOTE, VMBLST AND VMHC BOTH DECREASE OVER COURSE OF PREGNANCY - EXTRA "PROTECTION" FROM GOING BELOW ZERO IS NEEDED
VMHCPvv = RSW(VMHCPvv+VMHCN .LT. .50*VMHCN, -0.50*VMHCN, VMHCPvv)
VMHHCOvv = RSW(VMHHCOvv+VMHHON .LT. .01*VMHHON, -0.99*VMHHON, VMHHCOvv)
VMCBCPvv = RSW(VMCBCPvv+VMBCP .LT. .01*VMBCP, -0.99*VMBCP, VMCBCPvv)
VMCBCOvv = RSW(VMCBCOvv+VMBCO .LT. .01*VMBCO, -0.99*VMBCO, VMCBCOvv)
VINTTvv = RSW(VINTTvv+VINTT .LT. .01*VINTT, -0.99*VINTT, VINTTvv)
VINTOvv = RSW(VINTOvv+VINTO .LT. .01*VINTO, -0.99*VINTO, VINTOvv)

! set variable Vmax's
VMLSTN=VMLvv+VMLSTN
VMBLSTN=VMBLvv+VMBLSTN
VMHCN=VMHCPvv+VMHCN  
VMHHON=VMHHCOvv+VMHHON
VMBCP=VMCBCPvv+VMBCP
VMBCO=VMCBCOvv+VMBCO
VINTT=VINTTvv+VINTT
VINTO=VINTOvv+VINTO

!other variability based on point (LOCAL) senstivity analyses
CONSTANT KDRBCEv=0.0
CONSTANT KIRBCEv=0.0
CONSTANT HCTv=0.0
CONSTANT TRCRv=0.0
CONSTANT HECEv=0.0
CONSTANT KRHCRv=0.0
CONSTANT KSIv=0.0
CONSTANT KRRBCEv=0.0
! lognormal distributions
KDRBCEvv=exp(GAUSS(LOG(KDRBCE)-(SQRT(LOG((KDRBCEv)**2.+1))**2.)/2.,SQRT(LOG((KDRBCEv)**2.+1))))-KDRBCE
KIRBCEvv=exp(GAUSS(LOG(KIRBCE)-(SQRT(LOG((KIRBCEv)**2.+1))**2.)/2.,SQRT(LOG((KIRBCEv)**2.+1))))-KIRBCE
KRRBCEvv=exp(GAUSS(LOG(KRRBCE)-(SQRT(LOG((KRRBCEv)**2.+1))**2.)/2.,SQRT(LOG((KRRBCEv)**2.+1))))-KRRBCE
HCTvv=exp(GAUSS(LOG(HCTN)-(SQRT(LOG((HCTv)**2.+1))**2.)/2.,SQRT(LOG((HCTv)**2.+1))))-HCTN
TRCRvv=exp(GAUSS(LOG(TRCR)-(SQRT(LOG((TRCRv)**2.+1))**2.)/2.,SQRT(LOG((TRCRv)**2.+1))))-TRCR
KRHCRvv=exp(GAUSS(LOG(KRHCR)-(SQRT(LOG((KRHCRv)**2.+1))**2.)/2.,SQRT(LOG((KRHCRv)**2.+1))))-KRHCR
HECEvv=exp(GAUSS(LOG(HECE)-(SQRT(LOG((HECEv)**2.+1))**2.)/2.,SQRT(LOG((HECEv)**2.+1))))-HECE
KSIvv=exp(GAUSS(LOG(KSI)-(SQRT(LOG((KSIv)**2.+1))**2.)/2.,SQRT(LOG((KSIv)**2.+1))))-KSI
!IF value goes below 1%, set to 1% of nominal mean
KDRBCEvv = RSW(KDRBCEvv+KDRBCE .LT. .01*KDRBCE, -0.99*KDRBCE, KDRBCEvv)
KIRBCEvv = RSW(KIRBCEvv+KIRBCE .LT. .01*KIRBCE, -0.99*KIRBCE, KIRBCEvv)
KRRBCEvv = RSW(KRRBCEvv+KRRBCE .LT. .01*KRRBCE, -0.99*KRRBCE, KRRBCEvv)
HCTvv=RSW(HCTvv+HCTN .LT. .01*HCT, -0.99*HCTN, HCTvv)
TRCRvv = RSW(TRCRvv+TRCR .LT. .01*TRCR, -0.99*TRCR, TRCRvv)
KRHCRvv = RSW(KRHCRvv+KRHCR .LT. .01*KRHCR, -0.99*KRHCR, KRHCRvv)
HECEvv = RSW(HECEvv+HECE .LT. .01*HECE, -0.99*HECE, HECEvv)
KSIvv = RSW(KSIvv+KSI .LT. .01*KSI, -0.99*KSI, KSIvv)
! set variable params
KDRBCE=KDRBCEvv+KDRBCE
KIRBCE=KIRBCEvv+KIRBCE
KRRBCE=KRRBCEvv+KRRBCE
HCTN=HCTvv+HCTN
TRCR=TRCRvv+TRCR
HECE=HECEvv+HECE
KSI=KSIvv+KSI
KOSI=KSIvv+KOSI
KRHCR=KRHCRvv+KRHCR

! Blood flow to each compartment (L/h)
QF=(QFC+QFCv)*VF 
QH=(QHC+QHCVv)*VH
QR=(QRC+QRCv)*VR
QS=QSC*VS+QSCv 
QB=QBC*VB
QD=QDC*VD
QC=QF+QH+QR+QS+QB+QD+QPLA	! Total flow

! TIMING COMMANDS
CONSTANT TSTOP=1		! STOP SIMULATION (HR)


END		! END OF INITIAL

!######################################################################
DYNAMIC

!GROWING PREG COMPARTMENTS

IF(PREG.EQ.0)THEN
PAFC=0
ELSE
PAFC=0.1
END IF

IF(PREG.EQ.1)THEN
VUTR=(80+8.2931*GWEEK+0.3546*GWEEK^2)/1000  !Abduljalil et al 2012
ELSE
VUTR=1E-99
END IF

IF (PREG.EQ.1)THEN
VFI=VFN+((0.1305*GWEEK)+(0.0008*GWEEK^2)) !Abduljalil et al 2012
ELSE
VFI=0
END IF
IF (PREG.EQ.1)THEN
VFINIT=((VFN+0.1305*GWEEKSTART+0.0008*GWEEKSTART^2))
ELSE
VFINIT=0
END IF

IF (PREG.EQ.1)THEN
VBLI=VBLN-0.003741*GWEEK+0.004375*GWEEK^2-8.748e-5*GWEEK^3!+0.0098*GWEEK !Abduljalil et al 2012
ELSE
VBLI=0
END IF
IF (PREG.EQ.1)THEN
VBLINIT=VBLN-0.003741*GWEEKSTART+0.004375*GWEEKSTART^2-8.748e-5*GWEEKSTART^3!+0.0098*GWEEK
ELSE
VBLINIT=0
END IF


IF (PREG.EQ.1)THEN
VSI=VSN-0.2313*GWEEK+0.02311*GWEEK^2-0.0003376*GWEEK^3 !Abduljalil et al 2012 - 2/3 OF FFM
ELSE
VSI=0
END IF
IF (PREG.EQ.1)THEN
VSINIT=VSN-0.2313*GWEEKSTART+0.02311*GWEEKSTART^2-0.0003376*GWEEKSTART^3
ELSE
VSINIT=0
END IF

IF (PREG.EQ.1)THEN
VRI=VRN+0.004209*GWEEK+0.002854*GWEEK^2-7.610E-5*GWEEK^3!!+((80+8.2931*GWEEK+0.3546*GWEEK^2)/1000)!FIT TO TOTAL FFM - 1/3 - UTIRUS
ELSE
VRI=0
END IF
IF (PREG.EQ.1)THEN
VRINIT=VRN+0.004209*GWEEKSTART+0.002854*GWEEKSTART^2-7.610E-5*GWEEKSTART^3
ELSE
VRINIT=0
END IF



IF (PREG.EQ.0)THEN
PFC=PFCN
ELSE
PFC=PFCN-1.3*GWEEK-0.1249*GWEEK^2+0.00273*GWEEK^3
END IF

IF (PREG.EQ.0)THEN
PFO=PFON
ELSE
PFO=PFON-0.365*GWEEK-0.039*GWEEK^2+0.000849*GWEEK^3
END IF


IF (PREG.EQ.0)THEN
HCT=HCTN
ELSE
HCT=HCTN-((0.0544*GWEEK+0.0021*GWEEK^2)/100)
END IF

IF (PREG.EQ.0)THEN
VMHHO=VMHHON
ELSE
VMHHO=VMHHON+0.68*GWEEK+0.18*GWEEK^2
END IF

IF (PREG.EQ.0)THEN
VMHC=VMHCN
ELSE
VMHC=VMHCN-2.29*GWEEK-0.279*GWEEK^2
END IF

IF (PREG.EQ.0)THEN
VMLST=VMLSTN  !20% LOSS OVER PREGNANCY - CONSERVATIVE, Weitman et al found no hepatic difference in Peroxinase activity in liver
ELSE
VMLST=VMLSTN-346*GWEEK+0.9703*GWEEK^2
END IF
VMBLN=VMBL1/(1+exp((VMBL2-AGE)/VMBL3))		! Blood Vmaxc CPF-oxon -> TCPy (umol/hr/L plasma), Dynamic: logistic func MOVED FOR PREG COMMENT THIS OUT FOR BOOTSTRAP
IF (PREG.EQ.1)THEN
VMBLST=VMBLSTN-1945*GWEEK+35.04*GWEEK^2!+3.12*GWEEK^3
ELSE
VMBLST=VMBLSTN
END IF
!VMBLST=VMBL*VBL*(1-HCT)	 ! Blood CPF-oxon -> TCPy


! Age
IF (SPECIES .EQ. 0) THEN
AGE=AGE0+TIME/24.0		! Rat age (D)
ELSE
AGE=AGE0+TIME/8766		! Human age (Y)
END IF

! Age dependent B4 & B8 for humans (Luecke 2007)
IF ((AGE .GE. AGE1) .AND. (AGE .LT. AGE2) .AND. (SPECIES .EQ. 1)) THEN
B4x=B4;B8x=0
ELSE IF ((AGE .GE. AGE3) .AND. (SPECIES .EQ. 1)) THEN
B4x=0;B8x=B8
ELSE
B4x=0;B8x=0
END IF

! BWT (kg) Gompertz (Luecke 2007)
IF (AGE .GE. AGE2) THEN					
BWTINIT=B5*EXP((B6/B7)*(1-EXP(-B7*(AGE-AGE2))))+B8x*(AGE-AGE3)		
ELSE
BWTINIT=B1*EXP((B2/B3)*(1-EXP(-B3*AGE)))+B4x*(AGE-AGE1)
END IF
BWT=(BWTINIT*(1.-BWSW)+VFET+VPLA+VUTR+(VFI-VFINIT)+(VBLI-VBLINIT)+(VSI-VSINIT)+(VRI-VRINIT))+BWST*BWSW

BWTG=BWTINIT*1000		! BWT kg to grams for compartment vol calcs

IF (SPECIES .EQ. 0) THEN
LBWT=LOG(BWTINIT)		! Log of BWT for rat brain calc
ELSE
LBWT=BWTINIT*1000
END IF

IF (PREG.EQ.0)THEN
VF=VFN
ELSE
VF=VFI
END IF

IF (PREG.EQ.0)THEN
VBL=VBLN
ELSE
VBL=VBLI
END IF


IF (PREG.EQ.0)THEN
VS=VSN
ELSE
VS=VSI
END IF

IF (PREG.EQ.0)THEN
VR=VRN
ELSE
VR=VRI
END IF

TSA=71.81*(BWT**0.425)*(HT**0.725)  !DUBOIS AND DUBOIS 1916 HUMAN SA - MOVED TO INCLUDE INCREASE DURING PREGNANCY
SA=SAF*TSA

IF (PREG.EQ.0)THEN
QR=(QRC+QRCv)*VR
ELSE
QR=((QRC+0.006)+QRCv)*VR
END IF

IF (PREG.EQ.0) THEN
VFET=1e-99
ELSE  
 VFET = 3.55 * (EXP(-16.081*EXP(-5.67E-4*GTIME))+ EXP(-140.178*EXP(-7.01E-4*GTIME)))
 !VFET = 0.01*EXP(13.6*(1-EXP(-0.0702*GWEEK)))
END IF

IF (PREG.EQ.0) THEN
VPLA=1e-99
ELSE
  !VPLA = 0.85*EXP(-9.434*EXP(-5.23E-4*GTIME))
  VPLA = (1e-19+0.00675*GWEEK+0.000258*GWEEK^2)
END IF


! Human Blood flow to placenta (L/hr)
IF (PREG.EQ.0) THEN
QPLA=0
ELSE
   QPla = 58.5 * VPla
END IF     

IF (PREG.EQ.0) THEN
QUTR=0
ELSE
   QUTR = 62 * VPla
END IF     



! Scaled permeability-area product
       PAF = PAFC * (VFet**0.75)

! Compartment volumes (frac of BWT)
VBC=VB0+(VB1*LBWT)+(VB2*LBWT**2)+(VB3*LBWT**3)+(VB4*LBWT**4)+(VB5*LBWT**5)+(VB6*LBWT**6)
								! Brain
VHC=VH0+(VH1*BWTG)+(VH2*BWTG**2)+(VH3*BWTG**3)+(VH4*BWTG**4)	! Liver
VBLC=VBL0+(VBL1*BWTG)+(VBL2*BWTG**2)				! Blood
VFC=VF0+(VF1*BWTG)+(VF2*BWTG**2)+(VF3*BWTG**3)+(VF4*BWTG**4)+(VF5*LBWT**5)+(VF6*LBWT**6)	
								! Fat
VAC=VA0+(VA1*BWTG)+(VA2*BWTG**2)+(VA3*BWTG**3)+(VA4*BWTG**4)+(VA5*BWTG**5)+(VA6*BWTG**6)	
								! Adipose
VKC=VK0+(VK1*BWTG)+(VK2*BWTG**2)				! Kidney (in rapid)
VSPC=VSP0+(VSP1*BWTG)+(VSP2*BWTG**2) 				! Spleen (in rapid)
VLC=VL0+(VL1*BWTG)+(VL2*BWTG**2)				! Lung 
VGIC=VGI0+(VGI1*BWTG)+(VGI2*BWTG**2)+(VGI3*BWTG**3)		! Intestines (in rapid)
VMC=VM0+(VM1*BWTG)+(VM2*BWTG**2)+(VM3*BWTG**3)+(VM4*BWTG**4)			! Muscle (in slow)
VSKC=VSK0+(VSK1*BWTG)+(VSK2*BWTG**2)+(VSK3*BWTG**3)+(VSK4*BWTG**4)+(VSK5*BWTG**5)
								! Skin (in slow)
VFC=VFC+VFCv

VRC=VKC+VSPC+VLC+VGIC

VSC=VMC+VSKC+VBMC+(VAC-VFC)					! Slow

VSC=VSC+VSCu
VRC=VRC+VRCu
VBC=VBC+VBCu
VHC=VHC+VHCu
VFC=VFC+VFCu
VBLC=VBLC+VBLCu
VDC=VDC+VDCu
VINTC=VINTC+VINTCu

! Compartment volumes (L)
VB=VBC*BWT
VH=VHC*BWT
!VBL=VBLC*BWT
!VF=VFC*BWT
VRN=VRC*BWTINIT
!VS=VSC*BWT
VD=VDC*BWT
VINT=VINTC*BWT
VLU=VLC*BWT	!LUNG added 5.13

! Blood flow to each compartment (L/h)
QF=(QFC+QFCv)*VF 
QH=(QHC+QHCVv)*VH
QR=(QRC+QRCv)*VR
QS=QSC*VS+QSCv 
QB=QBC*VB
QD=QDC*VD
QC=QF+QH+QR+QS+QB+QD+QPLA	! Total flow

!mucociliator
FRLUGI=1-FRACIN		!FRACTIONAL transfer of CPF to GI tract - per O'Flaherty Chromium model
FRLUGIO=1-FRACIN   !FRACTIONAL transfer of OXON to GI tract
RLUGIO=FRLUGIO*CIO*QALV
 
! Age dependent metabolism (comment out for bootstrap)
!VMBL=VMBL1/(1+EXP((VMBL2-AGE)/VMBL3))		! Blood Vmaxc CPF-oxon -> TCPy (umol/hr/L plasma), dynamic logistic


! Scaling Vmaxc to Vmax (umol/h) 
VINTO=VINTOC*VINT		   	! Sm int. CPF -> CPF-oxon
VINTT=VINTTC*VINT			   ! Sm int. CPF -> TCPy
VINTOX=VINTOXC*VINT			 ! Sm int. CPF-oxon -> TCPy
VMBCO=VMCBCO*VB			     ! Brain CPF -> CPF-oxon
VMBCP=VMCBCP*VB			     ! Brain CPF -> TCPy
VMBCPT=VMBCPTC*VB
VMAXOLU=VMAXOLUC*VLU		! LUNG CPF->OXON
VMAXTLU=VMAXTLUC*VLU		! LUNG CPF->TCP
VMAXPLU=VMAXPLUC*VLU		!LUNG OXON->TCP



! Enzyme activity (umol/h)
! AChE
SBACH=BACHE*VB
SDACH=DACHE*VD
SHACH=HACHE*VH
SLUACH=LUACHE*VLU	
SBLACH=BLACHE*VBL*(1-HCT)
SRBACH=RBCHE*VBL*HCT
! CAE
SBRCE=BRCE*VB	
SDACE=DACE*VD		
SHECE=HECE*VH	
SLUECE=LUECE*VLU		
SPLOCE=PLOCE*VBL*(1-HCT)	
! BuChE
SBBUC=BBUCE*VB
SDBUC=DBUCE*VD	
SHBUC=HBUCE*VH
SHBUC=HBUCE*VH
SBLBUC=BLBUCE*VBL*(1-HCT)

! CALCULATION OF ESTERASE BINDING SITES (UMOLES) (MAXWELL, 1987)
IBCE=SBACH/TRCE		! TOTAL BRAIN ACHE (UMOLES)
IDCE=SDACH/TRCE		! TOTAL DIAPHRAGM ACHE (UMOLES)
IHCE=SHACH/TRCE		! TOTAL HEPATIC ACHE (UMOLES)
IBLCE=SBLACH/TRCE	! TOTAL PLASMA ACHE (UMOLES)
ILUCE=SLUACH/TRCE
IRBCE=SRBACH/TRCE	! Tot RBC AchE (umol)

IBBE=SBBUC/TRBE		! TOTAL BRAIN BUCHE (UMOLES)
IDBE=SDBUC/TRBE		! TOTAL DIAPHRAGM BUCHE (UMOLES)
IHBE=SHBUC/TRBE		! TOTAL HEPATIC BUCHE (UMOLES)
ILUBE=SLUBUC/TRBE
IBLBE=SBLBUC/TRBE	! TOTAL PLASMA BUCHE (UMOLES)

IBCR=SBRCE/TRCR		! TOTAL BRAIN CAE (UMOLES)
IDCR=SDACE/TRCR		! TOTAL DIAPHRAGM CAE (UMOLES)
IHCR=SHECE/TRCR		! TOTAL HEPATIC CAE (UMOLES)
ILUCR=SLUECE/TRCR
IBLCR=SPLOCE/TRCR	! TOTAL PLASMA CAE (UMOLES)


! CALCULATION OF ENZ SYN RATE (UMOLE/HR)
! LIVER
KSHCE=IHCE*KDHCE	! SYN RATE OF HEPATIC ACHE (UMOLE/HR)
KSHBE=IHBE*KDHBE	! SYN RATE OF HEPATIC BCHE (UMOLE/HR)
KSHCR=IHCR*KDHCR	! SYN RATE OF HEPATIC BCHE (UMOLE/HR)
! LUNG
KSLUCE=ILUCE*KDLUCE	! SYN RATE OF LUNG ACHE (UMOLE/HR)
KSLUBE=ILUBE*KDLUBE	! SYN RATE OF LUNG BCHE (UMOLE/HR)
KSLUCR=ILUCR*KDLUCR	! SYN RATE OF LUNG BCHE (UMOLE/HR)
! BRAIN
KSBCE=IBCE*KDBCE	! SYN RATE OF BRAIN ACHE (UMOLE/HR)
KSBBE=IBBE*KDBBE	! SYN RATE OF BRAIN BCHE (UMOLE/HR)
KSBCR=IBCR*KDBCR	! SYN RATE OF BRAIN BCHE (UMOLE/HR)
! DIAPHRAGM
KSDCE=IDCE*KDDCE	! SYN RATE OF DIAPHRAGM ACHE (UMOLE/HR)
KSDBE=IDBE*KDDBE	! SYN RATE OF DIAPHRAGM BCHE (UMOLE/HR)
KSDCR=IDCR*KDDCR	! SYN RATE OF DIAPHRAGM BCHE (UMOLE/HR)
! PLASMA
KSBLCE=IBLCE*KDBLCE	! SYN RATE OF PLASMA ACHE (UMOLE/HR)
KSBLBE=IBLBE*KDBLBE	! SYN RATE OF PLASMA BCHE (UMOLE/HR)
KSBLCR=IBLCR*KDBLCR	! SYN RATE OF PLASMA BCHE (UMOLE/HR)
! RBC
KSRBCE=KDRBCE*IRBCE	! Synthesis rate of RBC (umol/hr)

! Variable enzymology
VMBCP=VMCBCPvv+VMBCP
VMBCO=VMCBCOvv+VMBCO
VINTT=VINTTvv+VINTT
VINTO=VINTOvv+VINTO


DISCRETE DOSEX				! Oral dose code
ODOSE=ODOSE+DORAL(indx)*(BWT/MWC)*FA	! CPF dose molar conversion
SCDOSE=SCDOSE+DSC(indx)*(BWT/MWC)	! SC dose
ODOSEO=ODOSEO+DORALO(indx)*(BWT/MWO)*FA !OXON DOSE MOLAR CONVERSION
ADOSEO=ADOSEO+AMTODOSE(indx)*FA/MWO
ADOSEC=ADOSEC+AMTCDOSE(indx)*FA/MWC
indx=indx+1
SCHEDULE DOSEX .AT. DT(indx)
END



! Exposure Conditions Based on User Defined Initial Amounts of Chemical (mg)
     IF (concppm.EQ.0.0) THEN
     concmg=concmgm/1000						!CONvERT MG/M3 to mg/L
  ELSE
       CONCmg = CONCppm/24451.    !Convert ppm to mg/Liter!
  ENDIF
      
CONSTANT   DOSEINTERVAL=24					!TIME BETWEEN DAILY DOSES

! Exposure Conditions Based on User Defined Initial Amounts of Chemical (mg)
     IF (concppmO.EQ.0.0) THEN
     concmgO=concmgmO/1000						!CONvERT MG/M3 to mg/L
  ELSE
       CONCmgO = CONCppmO/24451.    !Convert ppm to mg/Liter!
  ENDIF
              
! Scaled Pulmonary Ventilation Rate 
CONSTANT QPC=0  !(L/hr) FOR 1 KG ANIMAL (Gentry et al)
    QALV= (QPC * (BWT**0.74))*0.7     
ALGORITHM IALG=2
TERMT(TIME .GE. TSTOP)	! STATEMENT TO STOP EXECUTION

!NOTE THE FOLLOWING CODE ADDED 4/24/14 - AFTER THE PUBLISHED MANUSCRIPT. iT DOES NOT AFFECT MODEL OUTPUT TO DATE.
VOLUME=0.654+(AGE*0.0663)+(AGE^2*-0.000664)  !this is for water dosing only, for use by EPA to dose over time. Water consumption changes with age, this is fit to 90%
!data from EPA exposures factors handbook (2011), table 3-7 using Graphpad Prism
CONSTANT EVENTS=1 !TIMES/DAY TO DRINK
DRINKING=VOLUME/EVENTS   !NOTE THIS IS FOR A DRINKING WATER SCENARIO ONLY, DATA IS INPUT USING A EXCEL TABLE AND THE NUMBER OF TIMES/DAY THAT A PERSON DRINKS IS USED TO ESTIMATE DOSE FROM WATER

!######################################################################
DERIVATIVE

!for pregnancy model - tsp 2.15
 Days = (TIME / 24.0) + GDstart

IF (PREG .EQ. 1) THEN	
Gtime = TIME + GDstart*24.0  !time of gestation (hr)
ELSE
GTIME=0
END IF
GWEEK = GTIME/168
GWEEKSTART=GDSTART/7

! IV DOSE/INFUSION 
IF ((TIME .LE. IVE) .AND. (TIME .GE. IVS)) THEN
IVDOSE=IVR
ELSE 
IVDOSE=0.0
END IF

! Inhalation Exposures
!pulse instructions
!FOR BOTH INALAITON AND DERMAL, A PULSE IS USED TO TURN ON/OFF EXPOSURES HOURLY, DAILY, AND OR YEARLY
! THE FIRST (0,YEAR,W2IN)  IS THE LONG-TERM ON/OFF  - THE UNITS ARE days!!! SO, YEAR=365 WILL TURN OFF THE OVERALL EXPOSURE AFTER 1 YEAR - NO MORE PULSES WILL OCCUR AND NO MORE DOSING
! THE SECOND IS FOR WEEKELY CONTROL. AGAIN, UNITS ARE DAYS. D3IN (D3) WILL SET HOW MANY DAYS/WEEK DOSING OCCURS. 1-7 WILL WORK. 
! the final (pulse(s2IN,24,p2IN)) is how many hours/day - 1-24 can be set. Note, if 24 is set, then dosing will be continuous until d3(in) turns it off.
CIM = CONCMG*DUMMYIN
DUMMYIN=pulse(0,YEAR*24,W2IN*24)*pulse(0,7*24,d3IN*24)*(pulse(s2IN,24,p2IN))
CI=CIM/MWC*1000	!DOSE IN UMOL



!Dermal pulse - 
DUMMY = pulse(0,YEAR*24,W2*24)*pulse(0,7*24,d3*24)*(pulse(s2,24,p2))  !this is a change from the published model, 2015
		!original published multi-route model: pulse(0,7*24,d3*24)*(pulse(s2,24,p2)) 





! CENTRAL SC COMPARTMENT
RASC=(KPS*APSC)-(KSC*ASC)-(KSP*ASC)	! CHANGE IN SC COMPARTMENT (uMOL/HR)
ASC=INTEG(RASC,DSCDC)+SCDOSE		! AMT IN SC COMPARTMENT (uMOL)
RASCA=KSC*ASC				! RATE OF ABSORPTION (uMOL/HR)

! PERIPHERAL SC COMPARTMENT (PSC)
RAPSC=(KSP*ASC)-(KPS*APSC)		! CHANGE IN PSC COMPARTMENT (uMOL/HR)
APSC=INTEG(RAPSC,0.0)			! AMT IN PSC COMPARTMENT (uMOL)


! FEED DOSE 
! SWITCHES FOR FEED OFF/ON
IF ((TIME .LE. FEEDT)) THEN
FEEDDOSE=FEDOSE1
ELSE 
FEEDDOSE=0.0
END IF

! DERMAL ABSORPTION FROM SKIN
!!! ADDED CJH 20JUL2015...
! Rate of change of amount on skin now has three terms:
!  1. Depletion from absorption into skin, via KPL
!  2. Depletion from dissipative loss, via KDISS (1/day)
!  3. Accumulation from daily exposure, via DERM (mg/kg/day)
!
! Compute rate at which CPF is being deposited on skin
! DDOSE is amt deposited each 24 hours, so convert to a 1/h rate - THIS IS NOT 24 HR, IT IS DOSE DEPOSITED/time of exposure/day(TSP)
VMDERMO=VMCDERMCO*VSKCC		     	! DERMAL CPF -> CPF-oxon
VMDERMC=VMCDERMCP*VSKCC						! DERMAL CPF -> TCPy
VMDERMEST=VMDERMESTC*VSKCC			! DERMAL ESTERASE OXON -> TCPy

DDOSE=((DERM*BWT)/MWC)*FAD !total amount uMol use when dose is ug/kg/day 
ASURF=INTEG(RASURF,DDOSE)

CSURF=ASURF/VLIQ

RASURF=-((((KPL*SA/1000)*CSURF))*DUMMY) - (KDIS/24.0)*ASURF

RADL=((((KPL*SA/1000)*CSURF))*DUMMY)-RADLL-RMETDT-RMETDO


ADL=INTEG(RADL,0.0)

RADLL=ADL*PSKL					!RATE OF ABSORPTION
ADLT=INTEG(RADLL,0.0)   !TOT ABSORBED

CADERM=ADL/VSKCC

RMETDT=CADERM*VMDERMC/(KMDERMCP+CADERM)!METAB TO TCPY
RMETDO=CADERM*VMDERMO/(KMDERMCO+CADERM)!METAB TO OXON
ADLO=INTEG(RMETDO,0.0)

RDERMTOTOX=RMETDO-(PSKO*ADERMOXTOT)-RDERMOX	! Rate of CPF-oxon change in skin
ADERMOXTOT=INTEG(RDERMTOTOX,0.0)		! Amt of CPF-oxon in skin
RDERMTO=PSKO*ADERMOXTOT		! Rate of CPF-oxon transfer to blood
ADLTO=INTEG(RDERMTO,0.0)  ! TOTAL OXON ABSORBED
CADLO=ADERMOXTOT/VSKCC			! Conc of CPF-oxon in skin
ADERMCO=INTEG(RDERMTO,0.0)

RDERMOX=CADLO*VMDERMEST/(KMDERMEST+CADLO)

ADLTCPYC=INTEG(RMETDT,0.0)
ADLTCPYO=INTEG(RMETDO,0.0)
TOTDERMTCPY=ADLTCPYC+ADLTCPYO

!######################################################################
! Stomach
RSTOM=FEEDDOSE-KSI*ASTOM-KAS*ASTOM+AMUC*KMUC !NOTE RLUGI IS FROM LUNG "ELEVATOR", STRUCTURE BASED ON O'Flaherty Chromium model
ASTOM=ODOSE+ADOSEC+ODOSEC+INTEG(RSTOM,0.0)	! Amt of CPF in stomach
RSTOML=KAS*ASTOM
ASTOML=INTEG(RSTOML,0.0)
RSTOMO=-KOSI*ASTOMO
ASTOMO=ODOSEO+ADOSEO+INTEG(RSTOMO,0.0)

RLUGI=FRLUGI*CI*QALV-AMUC*KMUC
AMUC=INTEG(RLUGI,0.0)

! Intestine
RINTC=(KSI*ASTOM)-MINTO-MINTTC-(KAI*AINTC)	! Rate of CPF change in intestine (umol/h)
AINTC=INTEG(RINTC,0.0)     			! Amt of CPF in intestine (umol)
RAO=KAI*AINTC            			! Rate of CPF transfer to liver 
CINTC=AINTC/VINT				! Conc of CPF in intestine (umol/L)
MINTO=VINTO*CINTC/(KMINTO+CINTC)		! Metab CPF to CPF-oxon
MINTTC=VINTT*CINTC/(KMINTC+CINTC)		! Metab CPF to TCPy

RINTO=MINTO-(KOS*AINTO)-MINTOX+(KOSI*ASTOMO)	! Rate of CPF-oxon change in intestine
AINTO=INTEG(RINTO,0.0)		! Amt of CPF-oxon in intestine w/ CPF-oxon oral dose
RLIVO=AINTO*KOS				! Rate of CPF-oxon transfer to liver
CINTO=AINTO/VINT			! Conc of CPF-oxon in intestine
MINTOX=VINTOX*CINTO/(KMINTOX+CINTO)	! Metab CPF-oxon to TCPy (PON1)

RINTP=MINTTC+MINTOX-RAP			! Rate of TCPy change in intestine
AINTP=ODOSEP+INTEG(RINTP,0.0)		! Amt of TCPy in intestine w/ TCPy oral dose
CINTP=AINTP/VINT			! Conc of TCPy in intestine
RAP=KIHP*AINTP				! Rate of TCPy transfer from intestine to liver



!#######################################################################
        
! MODEL 1 CPF

! LIVER (C) - MODEL 1
RAH=QH*(CAF-CVHF)-RHCO-RHCP+RAO+RSTOML		! RATE OF CHANGE IN LIVER
AH=INTEG(RAH,0.0)					! AMT CPF IN LIVER
CH=AH/VH						! CONC CPF IN LIVER
CVHF=CH/PHC						! FREE CPF IN BLD LEAVING LIVER
HAUC=INTEG(CH,0.0)					! AUC CPF IN LIVER
RHCO=VMHHO*CH/(KMHHCO+CH)				! Rate of CPF -> CPF-oxon
AMOX=INTEG(RHCO,0.0)		!AMOUNT METABOLIZED TO OXON - UMOL
RHCP=VMHC*CH/(KMHCP+CH)					! Rate of CPF -> TCPy
AMTCP=INTEG(RHCP,0.0)  !AMOUNT METABOLIZED TO TCPY - UMOL

RALINH = QALV * FRACIN * CI*N !RATE INTO LUNG
RALU=QC*(CAF-CLU)-(QALV*CLEX) - RAMLUOX-RAMLUT + RALINH! RATE OF CHANGE IN LUNG
ALU=INTEG(RALU,0.0)	! AMT CPF IN LUNG
CLU=ALU/VLU		! CONC CPF IN LUNG
CVLU=CLU/PLU		! FREE CONC CPF IN BLOOD LEAVING LUNG
FAULU=INTEG(CLU,0.0)	! AUC CPF IN LUNG
ALUINH=INTEG(RALINH,0.0)!UMOL ABSORBED
TALUINH=ALUINH
INMG=ALUINH*MWC/BWT/1000  !MG/KG ABSORBED


AMLUGI=INTEG(RLUGI,0.0)
RMUCOS=AMLUGI*KMUC
ALUNGDOSE=INTEG(RMUCOS,0.0)

!  AMOUNT EXHALED 
    CLEX = CVLU / PLAIR           ! CONCENTRATION, UMOL/L
     RAEX = QALV * CLEX
    AEX = INTEG(RAEX, 0.0)   ! AMOUNT, MG PER  
    AEXC = AEX * N		   ! AMOUNT, MG, FOR A GROUP OF RATS


!----------LUNG METABOLISM----
  RAMLUOX = VMAXOLU*CLU/(KMOLU+CLU)!rate to oxon
 AMLUOX = INTEG(RAMLUOX, 0.0) 
CLUOX=AMLUOX/PLUO
  RAMLUT = VMAXTLU*CLU/(KMTLU+CLU) !rate to tcpy
        AMLUOT=INTEG(RAMLUT,0.0)

! FAT (C) - MODEL 1
RAF=QF*(CAF-CVFF)	! RATE OF CHANGE IN FAT
AF=INTEG(RAF,0.0)	! AMT CPF IN FAT
CF=AF/VF		! CONC CPF IN FAT
CVFF=CF/PFC		! FREE CONC CPF IN BLOOD LEAVING FAT
FAUC=INTEG(CF,0.0)	! AUC CPF IN FAT


! SLOWLY PERF (C) - MODEL 1
RAS=QS*(CAF-CVSF)	! RATE OF CHANGE IN SLOW PERF
AS=INTEG(RAS,0.0)	! AMT CPF IN SLOW PERF
CS=AS/VS		! CONC CPF IN SLOW PERF
CVSF=CS/PSC		! FREE CPF IN BLOOD LEAVING SLOW PERF


! DIAPHRAM (C)- MODEL 1
RAD=QD*(CAF-CVDF)-RDPCO-RDPCP	! RATE OF CHANGE IN DIAPHRAM
AD=INTEG(RAD,0.0)	! AMT CPF IN DIAPHRAM
CD=AD/VD		! CONC CPF IN DIAPHRAM
CVDF=CD/PDC		! FREE CONC CPF IN BLOOD LEAVING DIAPHRAM
RDPCO=VMCDPO*CD/(KMDPCO+CD)	! Rate of CPF -> CPF-oxon
RDPCP=VMCDPP*CD/(KMDPCP+CD)	! Rate of CPF -> TCPy
ADCO=INTEG(RDPCO, 0.0)


! RAPID PERF (C) - MODEL 1
RAR=QR*(CAF-CVRF)	! RATE OF CHANGE IN RAPID PERF
AR=INTEG(RAR,0.0)	! AMT CPF IN RAPID PERF
CR=AR/VR		! CONC CPF IN RAPID PERF
CVRF=CR/PRC		! FREE CPF IN BLOOD LEAVING RAPID PERF


! BRAIN (C) 
RAB=QB*(CAF-CVBF)-RBCO-RBCP     ! RATE OF CHANGE IN BRAIN 
AB=INTEG(RAB,0.0)		! AMT CPF IN BRAIN 
CB=AB/VB                        ! CONC CPF IN BRAIN 
CVBF=CB/PBC                     ! FREE CONC CPF IN BLOOD LEAVING BRAIN 
RBCO=VMBCO*CB/(KMBCO+CB)	! Rate of CPF -> CPF-oxon
RBCP=VMBCP*CB/(KMBCP+CB)	! Rate of CPF -> TCPy

! Amount in Fetus (mg)
     RAFet = PAF * (CPla - CFet)
      AFet = INTEG(RAFet, 0.0)
      CFet = AFet / VFet
CFETPG=CFET*350.6*1000
   AUCCFet = INTEG(CFet, 0.0)
! Amount in Placenta (mg)
     RAPla = (QPla * (CAF - CVPla)) + (PAF * (CFET - CPLA))  !NOTE, ONE WAY TO GET MORE TANSFER WOULD BE TO USE TOTAL, NOT BOUND (CA INSTEAD OF CAF), BUT NEED TO DESCRIBE CORRECTLY
      APla = INTEG(RAPla, 0.0)
      CPla = APla / VPla
     CVPla = CPla / PPla
CVPLAPM=CVPla*350.6*1000


! BLOOD VENOUS & ARTERIAL (C) - MODEL 1
RABL=QC*(CV-CA)+RASCA+IVDOSE+RADLL	! RATE OF CHANGE IN MIXED BLOOD
RLBLD=QC*(CLU-CAF)  !LUNG BLOOD
ALBLD=INTEG(RLBLD,0.0)
ABL=INTEG(RABL,0.0)+INTEG(RLBLD,0.0)		! AMOUNT CPF IN MIXED BLOOD
CBL=ABL/VBL		! CONC CPF IN MIXED BLOOD
CBLB=CBL*FBC			! CONC CPF BOUND IN MIXED BLOOD
CBLF=CBL*(1-FBC)		! CONC FREE CPF IN ARTERIAL BLOOD
CV=CVEF+CBLB			! TOTAL CPF IN VENOUS BLOOD
CVPM=CV*350.6*1000
! TOTAL FREE CPF VENOUS BLOOD FROM TISSUES
CVEF=(QH*CVHF+QF*CVFF+QS*CVSF+QR*CVRF+QB*CVBF+QD*CVDF+CVPLA*QPLA)/QC
CA=CBL/PBL		! CONC CPF IN ARTERIAL BLOOD
CAPM=CA*350.6*1000 !conc cpf in arterial blood units of pg/ml
CAF=CA*FUBC !(1-FBC)			! CONC FREE CPF IN ARTERIAL BLOOD
BLAUC=INTEG(CBL,0.0)		! AUC CPF IN BLOOD
CAAUC=INTEG(CA,0.0)
BLAUCPM=INTEG(CAPM,0.0)
BLAUCAVG=BLAUC/TSTOP


!#######################################################################
! MODEL 2 CPF-OXON

! LIVER (O) - MODEL 2
RAHO=RHCO+QH*(CAOF-CVHOF)-HAE-(RHPCE+RHPBE+RHPCR)+RLIVO	! RATE OF CHANGE IN LIVER
AHO=INTEG(RAHO,0.0)					! AMT O IN LIVER
CHO=AHO/VH						! CONC O IN LIVER
CVHOF=CHO/PHO						! O IN BLD LEAVING LIVER
HAUCO=INTEG(CHO,0.0)					! AUC OXON IN LIVER

! HEPATIC A-ESTERASE ACTIVITY
HAE=VMLST*CHO/(KMLST+CHO)		! FORMATION OF TCP FROM OXON
AMPONH=INTEG(HAE,0.0)  !CLEARANCE OF OXN IN LIVER - UMOL
RHTCP=(RHPCE+RHPBE+RHPCR)+HAE+RHCP	! RATE OF TCP IN LIVER
AHTCP=INTEG(RHTCP,0.0)			! AMT TCP IN LIVER
CHTCP=AHTCP/VH				! CONC TCP IN LIVER
HAUCCP=INTEG(CHTCP,0.0)			! AUC TCP IN LIVER
                                
! HEPATIC B-ESTERASE ACTIVITY
HBE=AHCE+AHBE+AHCR			! TOTAL B-ESTERASE (UMOLES) = CHOLINE + BUTYRYL + CARBOXYL
IHE=100*(HBE/(ANHCE+ANHBE+ANHCR))	! % TOTAL B-ESTERASE INHIBITION
HCE=100*(AHCE/ANHCE)			! % CHOLINE-ESTERASE INHIBITION
HBES=100*(AHBE/ANHBE)			! % BUTURYL-ESTERSE INHIBITION
HCR=100*(AHCR/ANHCR)			! % CARBOXYL-ESTERASE INHIBITION

! CHOLINE ESTERASE ACTIVITY
RHCE=KSHCE-AHCE*(KDHCE+KIHCE*CHO)+HOCE*KRHCE	! RATE OF CHE ESTERASE
AHCE=INTEG(RHCE,IHCE)				! LIVER CHE ESTERASE (UMOLES)
RNHCE=KSHCE-ANHCE*KDHCE				! Rate of naive ChE 
ANHCE=INTEG(RNHCE,IHCE)				! Amt of naive ChE
RHOCE=AHCE*KIHCE*CHO-HOCE*(KAHCE+KRHCE)		! RATE OF CHE INHIBITION
HOCE=INTEG(RHOCE,0.0)				! LIVER CHE INHIBITED (UMOLES)
RHPCE=AHCE*KIHCE*CHO				! RATE OF TCP FORMATION
        
! BUTYRYL ESTERASE ACTIVITY
RHBE=KSHBE-AHBE*(KDHBE+KIHBE*CHO)+HOBE*KRHBE	! RATE OF BE ESTERASE ACT
AHBE=INTEG(RHBE,IHBE)				! LIVER BE ESTERASE ACT (UMOLES)
RNHBE=KSHBE-ANHBE*KDHBE				! Rate of naive BE act
ANHBE=INTEG(RNHBE,IHBE)				! Amt of naive BE
RHOBE=AHBE*KIHBE*CHO-HOBE*(KAHBE+KRHBE)		! RATE OF BE INHIBITION
HOBE=INTEG(RHOBE,0.0)				! LIVER BE INHIBITED (UMOLES)
RHPBE=AHBE*KIHBE*CHO				! RATE OF TCP FORMATION
        
! CARBOXYL ESTERASE ACTIVITY
RHCR=KSHCR-AHCR*(KDHCR+KIHCR*CHO)+HOCR*KRHCR	! RATE OF CR ESTERASE ACT
AHCR=INTEG(RHCR,IHCR)				! LIVER CR ESTERASE ACT (UMOLES)
RNHCR=KSHCR-ANHCR*KDHCR				! Rate of naive CR
ANHCR=INTEG(RNHCR,IHCR)				! Amt of naive CR
RHOCR=AHCR*KIHCR*CHO-HOCR*(KAHCR+KRHCR)		! RATE OF CR INHIBITION
HOCR=INTEG(RHOCR,0.0)				! LIVER CR INHIBITED (UMOLES)
RHPCR=AHCR*KIHCR*CHO				! RATE OF TCP FORMATION
      
!LUNG
CIMO = concmgO* pulse(0,7*24,5*24)*(pulse(s2,24,p2))
! for a 5 day/wk exposure, change first pulse to pulse(0,7*24,5*24)
! for daily, pulse(0,1e6,24)
CIO=CIMO/MWO
RALINHO = QALV * FRACIN * CIO*N !RATE INTO LUNG
RALUO=RAMLUOX+QALV*(CAOF-CVLUO) - LAE- (QALV*CLEX)+ RALINHO-(RLUPCE+RLUPBE+RLUPCR)! RATE OF CHANGE IN LUNG
ALUO=INTEG(RALUO,0.0)	! AMT CPF IN LUNG
CLUO=ALUO/VLU		! CONC CPF IN LUNG
CVLUO=CLUO/PLU		! FREE CONC CPF IN BLOOD LEAVING LUNG
FAULUO=INTEG(CLUO,0.0)	! AUC CPF IN LUNG
ALUINHO=INTEG(RALINHO,0.0)
TALUINHO=ALUINHO

!  AMOUNT EXHALED 
    CLEXO = CVLUO / PLAIRO          ! CONCENTRATION, MG/L
     RAEXO = QALV * CLEXO
    AEXO = INTEG(RAEXO, 0.0)   ! AMOUNT, MG PER  
    AEXCO = AEXO * N		   ! AMOUNT, MG, FOR A GROUP OF RATS


! PULMONARY A-ESTERASE ACTIVITY
LAE=VMAXPLU*CLUO/(KMPLUC+CLUO)		! FORMATION OF TCP FROM OXON
RLUTCP=(RLUPCE+RLUPBE+RLUPCR)+LAE+RAMLUT	! RATE OF TCP IN LUNG
ALUTCP=INTEG(RLUTCP,0.0)			! AMT TCP IN LUNG
CLUTCP=ALUTCP/VLU				! CONC TCP IN LUNG
LUAUCCP=INTEG(CLUTCP,0.0)			! AUC TCP IN LUNG
                                
! PULMONARY B-ESTERASE ACTIVITY
LUBE=ALUCE+ALUBE+ALUCR			! TOTAL B-ESTERASE (UMOLES) = CHOLINE + BUTYRYL + CARBOXYL
ILUE=100*(LUBE/(ANLUCE+ANLUBE+ANLUCR))	! % TOTAL B-ESTERASE INHIBITION
LUCE=100*(ALUCE/ANLUCE)			! % CHOLINE-ESTERASE INHIBITION
LUBES=100*(ALUBE/ANLUBE)			! % BUTURYL-ESTERSE INHIBITION
LUCR=100*(ALUCR/ANLUCR)			! % CARBOXYL-ESTERASE INHIBITION

! CHOLINE ESTERASE ACTIVITY
RLUCE=KSLUCE-ALUCE*(KDLUCE+KILUCE*CLUO)+LUOCE*KRLUCE	! RATE OF CHE ESTERASE
ALUCE=INTEG(RLUCE,ILUCE)				! LUNG CHE ESTERASE (UMOLES)
RNLUCE=KSLUCE-ANLUCE*KDLUCE				! Rate of naive ChE 
ANLUCE=INTEG(RNLUCE,ILUCE)				! Amt of naive ChE
RLUOCE=ALUCE*KILUCE*CLUO-LUOCE*(KALUCE+KRLUCE)		! RATE OF CHE INHIBITION
LUOCE=INTEG(RLUOCE,0.0)				! LUNG CHE INHIBITED (UMOLES)
RLUPCE=ALUCE*KILUCE*CLUO				! RATE OF TCP FORMATION
        
! BUTYRYL ESTERASE ACTIVITY
RLUBE=KSLUBE-ALUBE*(KDLUBE+KILUBE*CLUO)+LUOBE*KRLUBE	! RATE OF BE ESTERASE ACT
ALUBE=INTEG(RLUBE,ILUBE)				! LUNG BE ESTERASE ACT (UMOLES)
RNLUBE=KSLUBE-ANLUBE*KDLUBE				! Rate of naive BE act
ANLUBE=INTEG(RNLUBE,ILUBE)				! Amt of naive BE
RLUOBE=ALUBE*KILUBE*CLUO-LUOBE*(KALUBE+KRLUBE)		! RATE OF BE INHIBITION
LUOBE=INTEG(RLUOBE,0.0)				! LUNG BE INHIBITED (UMOLES)
RLUPBE=ALUBE*KILUBE*CLUO				! RATE OF TCP FORMATION
        
! CARBOXYL ESTERASE ACTIVITY
RLUCR=KSLUCR-ALUCR*(KDLUCR+KILUCR*CLUO)+LUOCR*KRLUCR	! RATE OF CR ESTERASE ACT
ALUCR=INTEG(RLUCR,ILUCR)				! LUNG CR ESTERASE ACT (UMOLES)
RNLUCR=KSLUCR-ANLUCR*KDLUCR				! Rate of naive CR
ANLUCR=INTEG(RNLUCR,ILUCR)				! Amt of naive CR
RLUOCR=ALUCR*KILUCR*CLUO-LUOCR*(KALUCR+KRLUCR)		! RATE OF CR INHIBITION
LUOCR=INTEG(RLUOCR,0.0)				! LUNG CR INHIBITED (UMOLES)
RLUPCR=ALUCR*KILUCR*CLUO				! RATE OF TCP FORMATION

! BRAIN OXON (O)- MODEL 2
RABO=QB*(CAOF-CVBOF)-(RBPCE+RBPBE+RBPCR)+RBCO-RBCOP	! RATE OF CHANGE IN BRAIN
ABO=INTEG(RABO,0.0)				! AMT O IN BRAIN
CBO=ABO/VB					! CONC O IN BRAIN
CVBOF=CBO/PBO					! CONC O IN BLD LEAVING BRAIN
BAUCO=INTEG(CBO,0.0)				! AUC O IN BRAIN
RBCOP=VMBCPT*CBO/(KMBCPT+CBO)	! Rate of CPF -> CPF-oxon


! BRAIN B-ESTERASE ACTIVITY
BBE=ABCE+ABBE+ABCR			! TOTAL B-ESTERASE = CHOLINE + BUTYRYL + CARBOXYL
IBE=100*(BBE/(ANBCE+ANBBE+ANBCR))	! % TOTAL B-ESTERASE INHIBITION
BCE=100*(ABCE/ANBCE)			! % CHOLINE-ESTERASE INHIBITION
BCEPERC=100-BCE						!PERCENT INHIBITED RATHER THAN INHIBITION
BBES=100*(ABBE/ANBBE)			! % BUTURYL-ESTERSE INHIBITION
BCR=100*(ABCR/ANBCR)			! % CARBOXYL-ESTERASE INHIBITION

! CHOLINE ESTERASE ACTIVITY
RBCE=KSBCE-ABCE*(KDBCE+KIBCE*CBO)+BOCE*KRBCE	! RATE OF CHE ESTERASE ACT
ABCE=INTEG(RBCE,IBCE)				! BRAIN CHE ESTERASE ACT (UMOLES)
RNBCE=KSBCE-ANBCE*KDBCE				! Rate of naive ChE act
ANBCE=INTEG(RNBCE,IBCE)				! Amt of naive ChE 
RBOCE=ABCE*KIBCE*CBO-BOCE*(KABCE+KRBCE)		! RATE OF CHE INHIBITION
BOCE=INTEG(RBOCE,0.0)				! BRAIN CHE INHIBITED (UMOLES)
RBPCE=ABCE*KIBCE*CBO				! RATE OF TCP FORMATION
        
! BUTYRYL ESTERASE ACTIVITY
RBBE=KSBBE-ABBE*(KDBBE+KIBBE*CBO)+BOBE*KRBBE	! RATE OF BE ESTERASE ACT
ABBE=INTEG(RBBE,IBBE)				! BRAIN BE ESTERASE ACT (UMOLES)
RNBBE=KSBBE-ANBBE*KDBBE				! Rate of naive BE act
ANBBE=INTEG(RNBBE,IBBE)				! Amt of naive BE
RBOBE=ABBE*KIBBE*CBO-BOBE*(KABBE+KRBBE)		! RATE OF BE INHIBITION
BOBE=INTEG(RBOBE,0.0)				! BRAIN BE INHIBITED (UMOLES)
RBPBE=ABBE*KIBBE*CBO 				! RATE OF TCP FORMATION
        
! CARBOXYL ESTERASE ACTIVITY
RBCR=KSBCR-ABCR*(KDBCR+KIBCR*CBO)+BOCR*KRBCR	! RATE OF CR ESTERASE ACT
ABCR=INTEG(RBCR,IBCR)				! BRAIN CR ESTERASE ACT (UMOLES)
RNBCR=KSBCR-ANBCR*KDBCR				! Rate of naive CRE act
ANBCR=INTEG(RNBCR,IBCR)				! Amt of naive CRE
RBOCR=ABCR*KIBCR*CBO-BOCR*(KABCR+KRBCR)		! RATE OF CR INHIBITION
BOCR=INTEG(RBOCR,0.0)				! BRAIN CR INHIBITED (UM L-1)
RBPCR=ABCR*KIBCR*CBO				! RATE OF TCP FORMATION


! DIAPHRAM OXON (O)- MODEL 2
RADO=QD*(CAOF-CVDOF)-(RDPCE+RDPBE+RDPCR)+RDPCO	! RATE OF CHANGE IN DIAPHRAM
ADO=INTEG(RADO, 0.0)				! AMT O IN DIAPHRAM
CDO=ADO/VD					! CONC O IN DIAPHRAM
CVDOF=CDO/PDO					! O IN BLD LEAVING DIAPHRAM
DAUCO=INTEG(CDO,0.0)				! AUC O IN DIAPHRAM

! DIAPHRAM B-ESTERASE ACTIVITY
DBE=ADCE+ADBE+ADCR			! TOTAL B-ESTERASE = CHOLINE + BUTYRYL + CARBOXYL
DBE1=ADCE+ADBE				! TOTAL CHOLINE + BUTURYL
IDE=100*(DBE/(ANDCE+ANDBE+ANDCR))	! % TOTAL B-ESTERASE INHIBITION
IDE1=100*(DBE1/(ANDCE+ANDBE))		! % TOTAL CHOLINE + BUTYRYL INHIBITION
DCE=100*(ADCE/ANDCE)			! % CHOLINE-ESTERASE INHIBITION
DBES=100*(ADBE/ANDBE)			! % BUTURYL-ESTERSE INHIBITION
DCR=100*(ADCR/ANDCR)			! % CARBOXYL-ESTERASE INHIBITION

! AChE ACTIVITY
RDCE=KSDCE-ADCE*(KDDCE+KIDCE*CDO)+DOCE*KRDCE	! RATE OF CHE ESTERASE ACT
ADCE=INTEG(RDCE,IDCE)				! DIAPHRAGM CHE ESTERASE ACT (UMOLES)
RNDCE=KSDCE-ANDCE*KDDCE				! Rate of naive ChE act
ANDCE=INTEG(RNDCE,IDCE)				! Amt of naive ChE
RDOCE=ADCE*KIDCE*CDO-DOCE*(KADCE+KRDCE)		! RATE OF CHE INHIBITION
DOCE=INTEG(RDOCE,0.0)				! DIAPHRAGM CHE INHIBITED (UMOLES)
RDPCE=ADCE*KIDCE*CDO				! RATE OF TCP FORMATION
        
! BUTYRYL ESTERASE ACTIVITY
RDBE=KSDBE-ADBE*(KDDBE+KIDBE*CDO)+DOBE*KRDBE	! RATE OF BE ESTERASE ACT
ADBE=INTEG(RDBE,IDBE)				! DIAPHRAGM BE ESTERASE ACT (UMOLES)
RNDBE=KSDBE-ANDBE*KDDBE				! Rate of naive BE act
ANDBE=INTEG(RNDBE,IDBE)				! Amt of naive BE
RDOBE=ADBE*KIDBE*CDO-DOBE*(KADBE+KRDBE)		! RATE OF BE INHIBITION
DOBE=INTEG(RDOBE,0.0)				! DIAPHRAGM BE INHIBITED (UMOLES)
RDPBE=ADBE*KIDBE*CDO				! RATE OF TCP FORMATION
        
! CARBOXYL ESTERASE ACTIVITY
RDCR=KSDCR-ADCR*(KDDCR+KIDCR*CDO)+DOCR*KRDCR	! RATE OF CR ESTERASE ACT
ADCR=INTEG(RDCR,IDCR)				! DIAPHRAGM CR ESTERASE ACT (UMOLES)
RNDCR=KSDCR-ANDCR*KDDCR				! Rate of naive CR
ANDCR=INTEG(RNDCR,IDCR)				! Amt of naive CR
RDOCR=ADCR*KIDCR*CDO-DOCR*(KADCR+KRDCR)		! RATE OF CR INHIBITION
DOCR=INTEG(RDOCR,0.0)				! DIAPHRAGM CR INHIBITED (UMOLES)
RDPCR=ADCR*KIDCR*CDO				! RATE OF TCP FORMATION


! FAT (O) - MODEL 2
RAFO=QF*(CAOF-CVFOF)		! RATE OF CHANGE IN FAT
AFO=INTEG(RAFO,0.0)		! AMT O IN FAT
CFO=AFO/VF			! CONC O IN FAT
CVFOF=CFO/PFO			! FREE O IN BLD LEAVING FAT
FAUCO=INTEG(CFO,0.0)		! AUC O IN FAT

! AMOUNT IN FETUS (MG)
   RAFETO = PAF * (CPLAO - CFETO)  !note PAF for oxon = for CPF, estimated
     AFETO = INTEG(RAFETO, 0.0)
     CFETO = AFETO / VFET
   AUCCFETO = INTEG(CFETO, 0.0)
! AMOUNT IN PLACENTA (MG)
     RAPLAO = (QPLA * (CAOF - CVPLAO)) + (PAF * (CFETO- CPLAO))
      APLAO = INTEG(RAPLAO, 0.0)
     CPLAO = APLAO / VPLA
     CVPLAO = CPLAO/ PPLAO

! SLOWLY PERF (O) - MODEL 2
RASO=QS*(CAOF-CVSOF)		! RATE OF CHANGE IN SLOW PERF
ASO=INTEG(RASO,0.0)		! AMT O IN SLOW PERF
CSO=ASO/VS			! CONC O IN SLOW PERF
CVSOF=CSO/PSO			! O IN BLD LEAVING SLOW PERF
SAUCO=INTEG(CSO,0.0)		! AUC O IN SLOW PERF


! RAPID PERF (O) - MODEL 2
RARO=QR*(CAOF-CVROF)		! RATE OF CHANGE IN RAPID PERF
ARO=INTEG(RARO,0.0)		! AMT OXON IN RAPID PERF
CRO=ARO/VR			! CONC OXON IN RAPID PERF
CVROF=CRO/PRO			! FREE OXON IN BLD LEAVING RAPID PERF
RAUCO=INTEG(CRO,0.0)		! AUC OXON IN RAPID PERF


! BLOOD (O) - MODEL 2
RABLO=QC*(CVO+CVLUO-CAO)-RBLOP+RDERMTO	! RATE OF CHANGE IN  MIXED BLOOD COMPT
ABLO=INTEG(RABLO,0.0)!+INTEG(RLBLDO,0.0)			! AMT O IN MIXED BLOOD COMPT
CBLO=ABLO/VBL				! CONC O IN  MIXED BLOOD COMPT
CAO=CBLO/PBLO				! CONC O IN ARTERIAL BLOOD
CAOF=ABLO*(1-FBO)/VBL			! CONC OF FREE O IN ARTERIAL BLOOD
CAOB=ABLO*(FBO)/VBL			! Conc CPF-oxon bound in art. blood (umol/L)
BLAUCO=INTEG(CBLO,0.0)			! AUC O IN  VENOUS BLOOD
BLAUCOAVG = BLAUCO/TSTOP		! Added by C Tan to calculate average cpf-oxon conc
BLAUCOF=INTEG(CAOF,0.0)
RBLOPA=VMBLST*CAOF/(KMBLST+CAOF)	! A-esterase CPF-oxon -> TCPy
AMPONBL=INTEG(RBLOPA,0.0)				!AMOUNT OXON CLEARED BY PLASMA PON1
RBLOPB=RBLPBE+RBLPCR+RBLPCE+RRBPCE	! Total B-esterase CPF-oxon -> TCPy
RBLOP=RBLOPA+RBLOPB			! Total CPF-oxon -> TCPy
CVO=(QH*CVHOF+QF*CVFOF+QS*CVSOF+QR*CVROF+QB*CVBOF+QD*CVDOF+CVPLAO*QPLA)/QC+CAOB
					! CONC CPF-oxon VENOUS BLOOD
NMCVO=CVO*1000		!IN nM TO PLACE ON OUTPUT SCALE

! BLOOD/PLASMA B-ESTERASE ACTIVITY
BLBE=ABLCE+ABLBE+ABLCR				! TOTAL B-ESTERASE = AChE + BuChE + CaE
BLBE1=ABLCE+ABLBE				! TOTAL AChE + BuChE
IBLE=100*(BLBE/(ANBLCE+ANBLBE+ANBLCR))		! % Total B-esterase inhibition in plasma
IBLE1=100*(BLBE1/(ANBLCE+ANBLBE))		! % Total AChE + BuChE inhibition in plasma
IBLE1PERC=100-IBLE1						!PERCENT INHIBITED RATHER THAN INHIBITION
BLCE=100*(ABLCE/ANBLCE)				! % AChE inhibition in plasma
BLBES=100*(ABLBE/ANBLBE)			! % BuChE inhibition in plasma
BLCR=100*(ABLCR/ANBLCR)				! % CaE inhibition in plasma
RBCCE=100*(ARBCE/ANRBCE)			! % AChE inhibition in RBC
RBCPERC=100-RBCCE							!CALCULATE % INHIBITED INSTEAD OF INHIBITION
RBCAVG=INTEG(RBCPERC,0.0)/TSTOP !AVG ACTIVITY (NOT INHIBITION)
TBLCE=100*((ABLCE+ARBCE)/(ANBLCE+ANRBCE))	! % AChE inhibition in blood
        
! Plasma AChE Activity
RBLCE=KSBLCE-ABLCE*(KDBLCE+KIBLCE*CBLO)+BLOCE*KRBLCE	! RATE OF CHE ESTERASE ACT
ABLCE=INTEG(RBLCE,IBLCE)      				! VENOUS CHE ESTERASE ACT (UMOLES)
RNBLCE=KSBLCE-ANBLCE*KDBLCE				! Rate of venous naive ChE
ANBLCE=INTEG(RNBLCE,IBLCE)				! Amt of venous naive ChE
RBLOCE=ABLCE*KIBLCE*CBLO-BLOCE*(KABLCE+KRBLCE)		! RATE OF CHE INHIBITION
BLOCE=INTEG(RBLOCE,0.0)					! VENOUS CHE INHIBITED (UMOLES)
RBLPCE=ABLCE*KIBLCE*CBLO				! RATE OF TCP FORMATION
      
! Plasma BuChE Activity
RBLBE=KSBLBE-ABLBE*(KDBLBE+KIBLBE*CBLO)+BLOBE*KRBLBE	! RATE OF BE ESTERASE ACT
ABLBE=INTEG(RBLBE,IBLBE)				! VENOUS BE ESTERASE ACT (UMOLES)
RNBLBE=KSBLBE-ANBLBE*KDBLBE				! Rate of venous naive BE
ANBLBE=INTEG(RNBLBE,IBLBE)				! Amt of venous naive BE
RBLOBE=ABLBE*KIBLBE*CBLO-BLOBE*(KABLBE+KRBLBE)		! RATE OF BE INHIBITION
BLOBE=INTEG(RBLOBE,0.0)					! VENOUS BE INHIBITED (UMOLES)
RBLPBE=ABLBE*KIBLBE*CBLO				! RATE OF TCP FORMATION
        
! Plasma CaE Activity
RBLCR=KSBLCR-ABLCR*(KDBLCR+KIBLCR*CBLO)+BLOCR*KRBLCR	! RATE OF CR ESTERASE ACT
ABLCR=INTEG(RBLCR,IBLCR)				! VENOUS CR ESTERASE ACT (UMOLES)
RNBLCR=KSBLCR-ANBLCR*KDBLCR				! Rate of venous naive CR
ANBLCR=INTEG(RNBLCR,IBLCR)				! Amt of venous naive CR
RBLOCR=ABLCR*KIBLCR*CBLO-BLOCR*(KABLCR+KRBLCR)		! RATE OF CR INHIBITION
BLOCR=INTEG(RBLOCR,0.0)					! VENOUS CR INHIBITED (UMOLES)
RBLPCR=ABLCR*KIBLCR*CBLO				! RATE OF TCP FORMATION

! RBC AChE Activity
RRBCCE=KSRBCE-ARBCE*(KDRBCE+KIRBCE*CBLO)+RBCOCE*KRRBCE	! RATE OF RBC CHE ESTERASE ACT
ARBCE=INTEG(RRBCCE,IRBCE)				! VENOUS RBC CHE ESTERASE ACT (UMOLES)
RNRBCCE=KSRBCE-ANRBCE*KDRBCE				! Rate of RBC naive ChE
ANRBCE=INTEG(RNRBCCE,IRBCE)				! Amt of RBC naive ChE
RRBOCE=ARBCE*KIRBCE*CBLO-RBCOCE*(KARBCE+KRRBCE)		! RATE OF CHE INHIBITION
RBCOCE=INTEG(RRBOCE,0.0)				! VENOUS CHE INHIBITED (UMOLES)
RRBPCE=ARBCE*KIRBCE*CBLO				! RATE OF TCP FORMATION


! TOTAL AUC FOR OXON IN ALL COMPARTMENTS
TAUCO=HAUCO+BAUCO+DAUCO+SAUCO+RAUCO+FAUCO+BLAUCO 

!#######################################################################
! TCPy
! TOTAL RATE OF TCP FORMATION
RP=RHCP+HAE+RHPCE+RHPBE+RHPCR+RBPCE+RBPBE+RBPCR+RDPCE+RDPBE+RDPCR+RBLOP+RAP+RBCP+RLUTCP+RBCOP+RDPCP+RMETDT
ATCP=INTEG(RP,0.0)		! TCP (UMOLES)
RA1=RP-RAE			! Rate of change in TCPy compartment 1 (umol/h)
A1=INTEG(RA1,0.0)		! Amt of TCPy in compartment 1 (umol)
CBLP=A1/V1 			! TCPy conc in compartment 1 (umol/L)
RAE=A1*KE			! Rate of TCPy excretion (umol/h)
AE=INTEG(RAE,0.0)		! Amt of TCPy excreted (umol)
AEUG=AE*198.4
urinetcpy = RAE*198.4/VVOL           !conc in (ug/L), C.Tan
urinetumol=RAE/VVOL  !!conc in (umol/L), assuming urine flow rate = 2.4 L/h C. Tan
!Blanker et al, 2001 = 1.5 L/day  = 0.06 L/hr
!Hays et al., 2015 - 0.7 ml/hr/kg = 0.05 L/hr - LOW 0.02, HIGH 0.4


END     	! END OF DERIVATIVE

!#######################################################################
! Total mass of CPF, CPF-oxon, and TCPy (umol)              
TMASSC=ASTOM+AINTC+AH+AF+AS+AD+AR+AB+ABL+ASC+APSC+ALU+APLA+AFET	! CPF
TMASSO=AINTO+AHO+AFO+ASO+ADO+ARO+ABO+ABLO+ALUO+ADLO+APLAO	      ! CPF-oxon
TMASSP=AINTP+A1+AE				                       ! TCPy
TMASS=TMASSC+TMASSO+TMASSP			                 ! Total mass
! Mass Balance
TDOSE=ODOSE+ODOSEC+ADOSEC+ADOSEO+ODOSEP+ODOSEO+SCDOSE+TIV+DSCDC+ADLT+ALUINH+ALUNGDOSE-AEX-AEXO+FEEDBALANCE
TDOSEUG=TDOSE*MWC
MASSBAL=TDOSE/(TMASS+1E-9)					! Should ~1
!Tissue volume checks
VTOTCOMPARTMENTS=(VHC+VFC+VSC+VRC+VDC+VBC+VBLC+VLC)/0.9 !should ~0.9-1.1 
BWTotal=(0.9*BWT)-(VH+VF+VS+VD+VR+VB+VBL+VPLA+VUTR) !should <10 for a human

END     	! END OF DYNAMIC

END     	! END OF PROGRAM
