The Binding of Symmetric Ligands of Two Types to an Isotropic Linear Lattice
by Alan R. Wolfe

This is the initialization notebook for the isotropic lattice, 2 types of symmetic ligand case.  It needs to be evaluated only once per Mathematica  session.  After opening this notebook, you should get a dialog box asking whether you want to automatically evaluate all the initialization cells; click Yes (otherwise, execute the Action-Evaluate Initialization command).  The need to respond to the dialog box can be avoided by selecting the "Always" option for "On opening a Notebook, load initialization cells:" at the bottom of the Action Preferences box under the Edit menu.  After evaluation, close this notebook (without quitting Mathematica ) and open the output notebook.  You should  quit Mathematica  and start a new session when switching to a different case.

Turn off General::spell1 warning

Off [ General::spell1 ] ;

Establish Mathematica version number

versionIsLessThan2 = True ; 
If [ $VersionNumber >= 2, versionIsLessThan2 = False ] ;

Choose whether output is to be saved

saveOutput = False ;

Definitions

points::usage = "points is the number of positive r1 values for which
parameter values will be calculated.  It also equals the number of intervals
on the r1-axis of the plots." ; 
kL1::usage = "kL1 is the binding constant for ligand 1-lattice
interactions." ; 
kL2::usage = "kL2 is the binding constant for ligand 2-lattice
interactions." ; 
nL1::usage = "nL1 is the number of lattice residues occupied per bound
ligand 1." ; 
nL2::usage = "nL2 is the number of lattice residues occupied per bound
ligand 2." ; 
wL1::usage = "wL1 is the cooperativity parameter for ligand 1-ligand 1
interactions on an isotropic lattice." ; 
wL2::usage = "wL2 is the cooperativity parameter for ligand 2-ligand 2
interactions on an isotropic lattice." ; 
wL12::usage = "wL12 is the cooperativity parameter for ligand 1-ligand 2
(or vice versa) interactions on an isotropic lattice." ; 
r1::usage = "r1 is the ratio of bound ligand 1 to lattice residues (plotted
on the abscissa of a Scatchard plot)." ; 
r2::usage = "r2 is the ratio of bound ligand 2 to lattice residues (plotted
on the abscissa of a Scatchard plot)." ; 
R2::usage = "R2 is the ratio of bound ligand 2 to lattice residues (used
only in equations in which it is being solved for)." ; 
f::usage = "f is the fraction of lattice residues that are free." ; 
F::usage = "F is the fraction of lattice residues that are free
(used in equations being solved for R2)." ; 
x1::usage = "x1 is r1/f." ; 
x2::usage = "x2 is r2/f." ; 
r1L1f::usage = "r1L1f is r1 divided by the free concentration of ligand 1
(plotted on the ordinate of a Scatchard plot)." ; 
r2L2f::usage = "r2L2f is r2 divided by the free concentration of ligand 2
(plotted on the ordinate of a Scatchard plot)." ; 
EL1::usage = "EL1 is the neighbor-effect parameter for a ligand 1
(symmetric) binding site on an isotropic lattice." ; 
EL2::usage = "EL2 is the neighbor-effect parameter for a ligand 2
(symmetric) binding site on an isotropic lattice." ; 
r1Sat::usage = "r1Sat is the limiting value of r1 for a given value of r2,
i.e., the ratio of ligand 1 to lattice residues when the lattice is saturated." ; 
r2Sat::usage = "r2Sat is the limiting value of r2, i.e., the ratio of 
ligand 2 to lattice residues when the lattice is saturated (and no other
ligand is present)." ; 
el2::usage = "el2 is the neighbor-effect parameter for ligand 2
(used only in calculating r2 on a lattice without any bound ligand 1)." ; 
fracRSat2::usage = "fracRSat2 is the fraction of the maximum possible
ratio of bound ligand 2 to lattice residues that is present on the lattice." ; 
freeL2::usage = "freeL2 is the concentration of free ligand 2." ; 
totalL2::usage = "totalL2 is the total concentration of ligand 2." ; 
polymer::usage = "polymer is the total lattice residue concentration." ;

r1L1f := Re [ kL1 f  ff^(nL1-1) EL1^2 ] ; 
r2L2f := Re [ kL2 f  ff^(nL2-1) EL2^2 ] ;

Determination of the maximum value of r1

r1SatDetermine::usage = "A routine to determine the ratio of bound
ligand 1 to lattice residues at saturation." ; 
r1SatDetermine := { 
     If [ type==1, 
          If [ wL1 >0 && wL2 >0 && wL12 >0, r1Sat = (1-r2 nL2)/nL1 ] ; 
          If [ wL1==0 && wL2 >0 && wL12 >0, 
               r1Sat = Min [ (1-r2 nL2)/nL1, (1 - r2 (nL2-1))/(nL1+1) ] ] ; 
          If [ wL1 >0 && wL2==0 && wL12 >0, 
               If [ nL1==1, r1Sat = 1-r2 nL2, 
                    r1Sat = Min [ (1-r2 nL2)/nL1, 
                    (1 - r2 (nL2+1))/(nL1-1) ] ] ] ; 
          If [ wL1 >0 && wL2 >0 && wL12==0, 
               r1Sat = (1-r2 nL2)/nL1 ] ; 
          If [ wL1 >0 && wL2==0 && wL12==0, 
               r1Sat = (1-r2 (nL2+1) )/nL1 ] ; 
          If [ wL1==0 && wL2 >0 && wL12==0, 
               r1Sat = (1-r2 nL2)/(nL1+1) ] ; 
          If [ wL1==0 && wL2==0 && wL12 >0, 
               If [ nL1==1, r1Sat = (1-r2 (nL2-1))/2, 
                    r1Sat = Min [ (1-r2 (nL2-1))/(nL1+1), 
                    (1 - r2 (nL2+1))/(nL1-1) ] ] ] ; 
          If [ wL1==0 && wL2==0 && wL12==0, 
               r1Sat = (1-r2 (nL2+1) )/(nL1+1) ] ; 
          If [ r1Sat==(1-r2 nL2)/nL1, 
               fMin = 0, fMin = 1 - nL1 r1Sat - nL2 r2 ] ] ; 
     If [ type==2 || type==3, 
          If [ wL1==0, r1Sat = 1/(nL1+1), r1Sat = 1/nL1 ] ] } ;

Equations for calculating values of neighbor-effect parameters on an unsaturated lattice

neChoose::usage = "A routine to determine which set of equations
for neighbor-effect parameters should be used when r1 < r1Sat." ; 
neChoose := { 
     If [ type==1, If [ wL1==wL2==wL12, 
          ne := ne2 ; neChoice = "ne2", ne := ne1 ; neChoice = "ne1" ] ] ; 
     If [ type==2 || type==3, If [ wL1==wL2==wL12, 
          ne := ne4 ; neChoice = "ne4", ne := ne3 ; neChoice = "ne3" ] ] ; 
     If [ neChoice=="ne2", needCheck = False, needCheck = True ] } ;

setLimits::usage = "setLimits is a routine that determines the
largest and smallest physically possible values for the neighbor-effect
parameters." ; 
setLimits := { 
     neTol = 0.0000001 ; 
     r2min =  - neTol ; 
     r2max = r2Sat + neTol ; 
     EL1min = Min [ 1, wL1, wL12 ] - neTol ; 
     EL1max = Max [ 1, wL1, wL12 ] + neTol ; 
     EL2min = Min [ 1, wL2, wL12 ] - neTol ; 
     EL2max = Max [ 1, wL2, wL12 ] + neTol ; 
     cpTol = 0.000001 ; 
     cpMin = - cpTol ; 
     cpMax = 1 + cpTol } ;

checkLimits::usage = "checkLimits is a routine that checks whether
the calculated values of the neighbor-effect parameters and conditional
probabilities are within their permitted ranges, and sets flag = True
(indicating that a computational error has occurred) if one of them is outside
its permitted range." ; 
checkLimits := { 
     If [ flag==False, 
          If [ Head [ EL1] == Complex || Head [ EL2] == Complex || 
               Head [ r2] == Complex, 
               {EL1,EL2,r2} = Map [ Re, {EL1,EL2,r2} ] ; 
               r1Error = r1 ; flag = True ]  ] ; 
     If [ flag == False, 
          If [ r2<r2min || r2>r2max || EL1<EL1min || EL1>EL1max || 
               EL2<EL2min || EL2>EL2max || f<0 || f>1, 
               r1Error = r1 ; flag = True ]  ] ; 
     If [ flag==False && r1 < r1Sat, 
          fF = ff ; bL1fTest = fF/EL1 ; bL2fTest = fF/EL2 ; 
          fbL1Test = x1 fF/EL1 ; fbL2Test = x2 fF/EL2 ; 
          bL1bL1Test = x1 wL1 fF/EL1^2 ; 
          bL2bL2Test = x2 wL2 fF/EL2^2 ; 
          bL1bL2Test = x2 wL12 fF/(EL1 EL2) ; 
          bL2bL1Test = x1 wL12 fF/(EL1 EL2) ] ; 
     If [ flag==False && r1==r1Sat, 
          fF = 0 ; bL1fTest = BL1F ; bL2fTest = BL2F ; 
          fbL1Test = FBL1 ; fbL2Test = FBL2 ; 
          bL1bL1Test = BL1BL1 ; bL2bL2Test = BL2BL2 ; 
          bL1bL2Test = BL1BL2 ; bL2bL1Test = BL2BL1 ] ; 
     If [ flag==False,           
          If [ fF < cpMin || fF > cpMax || 
               bL1fTest < cpMin || bL1fTest > cpMax || 
               bL2fTest < cpMin || bL2fTest > cpMax || 
               fbL1Test < cpMin || fbL1Test > cpMax || 
               fbL2Test < cpMin || fbL2Test > cpMax || 
               bL1bL1Test < cpMin || bL1bL1Test > cpMax || 
               bL2bL2Test < cpMin || bL2bL2Test > cpMax || 
               bL1bL2Test < cpMin || bL1bL2Test > cpMax || 
               bL2bL1Test < cpMin || bL2bL1Test > cpMax, 
               r1Error = r1 ; flag = True ]  ] } ;

printWarning1::usage = "printWarning1 warns that a computational error
has been detected (used when flag = True)." ; 
printWarning1 := { 
Print [ "WARNING:  A COMPUTATIONAL ERROR WAS DETECTED AT r1 = ", r1Error//N, "." ] ; 
Print [ "One or more of the calculated neighbor-effect parameter or" ] ; 
Print [ "conditional probability values lies outside its permitted range." ] ; 
Print [ "Refer to the 'Read Me' file for possible remedies." ] } ;

printWarning2::usage = "printWarning2 warns that a computational error
has been detected (used when flag = True)." ; 
printWarning2 := { 
Print [ "WARNING:  A COMPUTATIONAL ERROR WAS DETECTED AT r1 = ", r1Error//N, "." ] ; 
Print [ "Refer to the first warning message for more info." ] } ;

printWarning3::usage = "printWarning3 warns that a computational error
has been detected at r1 = zero (used when r1ZeroError = True)." ; 
printWarning3 := { 
Print [ "WARNING:  A COMPUTATIONAL ERROR WAS DETECTED AT r1 = 0." ] ; 
Print [ "Try increasing the value of the 'loopNumber' parameter (in the closed" ] ; 
Print [ "cell just above), or decreasing the value of 'freeL2' (for type 2" ] ; 
Print [ "isotherms) or 'totalL2' (for type 3 isotherms)." ] } ;

ne1::usage = "ne1 gives equations for the neighbor-effect parameters
for an unsaturated lattice (used for type 1 isotherms when wL1, wL2
and wL12 are not the same)." ; 
ne1 := { 
     {EL1,EL2} = {eL1,eL2}/.FindRoot [ { 
          eL1 (1 + x1/eL1 + x2/eL2) == 1 + wL1 x1/eL1 + wL12 x2/eL2 , 
          eL2 (1 + x1/eL1 + x2/eL2) == 1 + wL2 x2/eL2 + wL12 x1/eL1 } , 
          {eL1,EL1}, {eL2,EL2}, AccuracyGoal ->7, 
          MaxIterations ->200, WorkingPrecision ->18 ]  } ;

ne1::usage = "ne1 gives equations for the neighbor-effect parameters
for an unsaturated lattice (used for type 1 isotherms when wL1, wL2
and wL12 are not the same)." ; 
(* Alternate Version *) 
ne1 := { 
     {EL1,EL2} = {eL1,eL2}/.FindRoot [ { 
          eL1 == (1 + wL1 x1/eL1 + wL12 x2/eL2)/(1 + x1/eL1 + x2/eL2) , 
          eL2 == (1 + wL2 x2/eL2 + wL12 x1/eL1)/(1 + x1/eL1 + x2/eL2) } , 
          {eL1,EL1}, {eL2,EL2}, AccuracyGoal ->7, 
          MaxIterations ->200, WorkingPrecision ->18 ]  } ;

ne2::usage = "ne2 gives the value of the neighbor-effect parameters
below saturation for type 1 isotherms when wL1 = wL2 = wL12." ; 
ne2 := EL1 = EL2 = ( 1-x1-x2 + Sqrt [ (1-x1-x2)^2 + 4 (x1+x2) wL1 ] )/2 ;

ne3::usage = "ne3 gives equations for the neighbor-effect parameters
for an unsaturated lattice (used for type 2 or 3 isotherms when wL1, wL2
and wL12 are not the same)." ; 
ne3 := { 
     {EL1,EL2,r2} = {eL1,eL2,R2}/.FindRoot [ { 
          eL1 (F + r1/eL1 + R2/eL2) == F + wL1 r1/eL1 + wL12 R2/eL2 , 
          eL2 (F + r1/eL1 + R2/eL2) == F + wL2 R2/eL2 + wL12 r1/eL1 , 
          R2 eL2^(-2) == FreeL2 kL2 F (F/(F + r1/eL1 + R2/eL2))^(nL2-1) , 
          F == 1 - nL1 r1 - nL2 R2 } , 
          {eL1,EL1}, {eL2,EL2}, {R2,r2}, {F,f}, AccuracyGoal ->7, 
          MaxIterations ->200, WorkingPrecision ->18 ] ; 
     f = 1 - nL1 r1 - nL2 r2 ; 
     x1 = r1/f ; 
     x2 = r2/f } ;

ne3::usage = "ne3 gives equations for the neighbor-effect parameters
for an unsaturated lattice (used for type 2 or 3 isotherms when wL1, wL2
and wL12 are not the same)." ; 
(* Alternate Version *) 
ne3 := { 
     {EL1,EL2,r2} = {eL1,eL2,R2}/.FindRoot [ { 
          eL1 == (F + wL1 r1/eL1 + wL12 R2/eL2)/(F + r1/eL1 + R2/eL2) , 
          eL2 == (F + wL2 R2/eL2 + wL12 r1/eL1)/(F + r1/eL1 + R2/eL2) , 
          R2 eL2^(-2) == FreeL2 kL2 F (F/(F + r1/eL1 + R2/eL2))^(nL2-1) , 
          F == 1 - nL1 r1 - nL2 R2 } , 
          {eL1,EL1}, {eL2,EL2}, {R2,r2}, {F,f}, AccuracyGoal ->7, 
          MaxIterations ->200, WorkingPrecision ->18 ] ; 
     f = 1 - nL1 r1 - nL2 r2 ; 
     x1 = r1/f ; 
     x2 = r2/f } ;

ne4::usage = "ne4 gives the value of the neighbor-effect parameters
below saturation for type 2 or 3 isotherms when wL1 = wL2 = wL12." ; 
ne4 := { 
     {EL1,EL2,r2} = {eL2,eL2,R2}/.FindRoot [ { 
          eL2 F + r1 + R2 == (eL2 F + wL1 r1 + wL1 R2)/eL2 , 
          R2 eL2^(-2) == FreeL2 kL2 F (F eL2/(F eL2 + r1 + R2))^(nL2-1), 
          F == 1 - nL1 r1 - nL2 R2 }, 
          {eL2,EL2}, {R2,r2}, {F,f}, AccuracyGoal ->7, 
          MaxIterations ->200, WorkingPrecision ->18 ] ; 
     f = 1 - nL1 r1 - nL2 r2 ; 
     x1 = r1/f ; 
     x2 = r2/f } ;

neL2a::usage = "neL2a gives the values of EL2 and r2 at a known
concentration of free ligand 2 (used only when r1 = 0 with type 2
isotherms)." ; 
neL2a := { 
     {EL2,X2} = {eL2,z}/.FindRoot [ { 
          eL2 (eL2 + z) == eL2 + z wL2, 
          z eL2^(-2) == FreeL2 kL2 (1 + z/eL2)^(1-nL2) }, 
          {eL2, EL2}, {z, X2}, AccuracyGoal->7, 
          MaxIterations ->200, WorkingPrecision ->18 ] ; 
     r2 = X2/(1+X2 nL2 )  } ;

neL2b::usage = "neL2b gives the values of EL2 and r2 at a known
concentration of free ligand 2 (used only when r1 = 0 with type 3
isotherms)." ; 
neL2b := { 
     {EL2,X2,r2} = {eL2,z,R2}/.FindRoot [ { 
          eL2 (eL2 + z) == eL2 + z wL2, 
          z eL2^(-2) == FreeL2 kL2 (1 + z/eL2)^(1-nL2), 
          z == R2/(1-nL2 R2) }, 
          {eL2, EL2}, {z, X2}, {R2, r2}, AccuracyGoal->7, 
          MaxIterations ->200, WorkingPrecision ->18 ]  } ;

startPrint::usage = "startPrint starts the printout of the Scatchard plot." ; 
startPrint := { 
     Print [ "r1         r1/L1f" ] ; 
     If [ r1ZeroError==False, 
          Print [ "0          ", r1L1finitial//N ] ] } ;

Expressions for conditional probabilities when r1 is below its maximum value

ff::usage = "ff is the probability that a randomly chosen free
lattice residue is bordered on a given side by another free residue." ; 
bL1f::usage = "bL1f is the probability that a randomly chosen
bL1 bound ligand end is bordered by a free lattice residue." ; 
bL2f::usage = "bL2f is the probability that a randomly chosen
bL2 bound ligand end is bordered by a free lattice residue." ; 
fbL1::usage = "fbL1 is the probability that a randomly chosen free
lattice residue is bordered on a given side by a bL1 end of a bound
ligand." ; 
fbL2::usage = "fbL2 is the probability that a randomly chosen
free lattice residue is bordered on a given side by a bL2 end of a bound
ligand." ; 
bL1bL1::usage = "bL1bL1 is the probability that a randomly
chosen bL1 end of a bound ligand is bordered by another bL1 bound
ligand end." ; 
bL2bL2::usage = "bL2bL2 is the probability that a randomly chosen
bL2 end of a bound ligand is bordered by another bL2 bound ligand end." ; 
bL1bL2::usage = "bL1bL2 is the probability that a randomly chosen
bL1 end of a bound ligand is bordered by a bL2 bound ligand end." ; 
bL2bL1::usage = "bL2bL1 is the probability that a randomly chosen
bL2 end of a bound ligand is bordered by a bL1 bound ligand end." ;

ff := 1 / (1 + x1/EL1 + x2/EL2) ; 
bL1f [ h_] := h[[8]]/ h[[2]] /;  h[[1]] > 0 ; (* = ff/EL1 *) 
bL1f [ h_] := BL1F /;  h[[1]]==0 ; 
bL2f [ h_] := h[[8]]/ h[[3]] /;  h[[1]] > 0 ; (* = ff/EL2 *) 
bL2f [ h_] := BL2F /;  h[[1]]==0 ; 
fbL1 [ h_] := h[[6]] h[[8]]/ h[[2]] /;  h[[1]] > 0 ; (* = x1 ff/EL1 *) 
fbL1 [ h_] := FBL1 /;  h[[1]]==0 ; 
fbL2 [ h_] := h[[7]] h[[8]]/ h[[3]] /;  h[[1]] > 0 ; (* = x2 ff/EL2 *) 
fbL2 [ h_] := FBL2 /;  h[[1]]==0 ; 
bL1bL1 [ h_] := h[[6]] wL1 h[[8]]/ h[[2]]^2 /;  h[[1]] > 0 ; 
     (* = x1 wL1 ff/EL1^2 *) 
bL1bL1 [ h_] := BL1BL1 /;  h[[1]]==0 ; 
bL2bL2 [ h_] := h[[7]] wL2 h[[8]]/ h[[3]]^2 /;  h[[1]] > 0 ; 
     (* = x2 wL2 ff/EL2^2 *) 
bL2bL2 [ h_] := BL2BL2 /;  h[[1]]==0 ; 
bL1bL2 [ h_] := h[[7]] wL12 h[[8]]/(h[[2]] h[[3]]) /;  h[[1]] > 0 ; 
     (* = x2 wL12 ff/(EL1 EL2) *) 
bL1bL2 [ h_] := BL1BL2 /;  h[[1]]==0 ; 
bL2bL1 [ h_] := h[[6]] wL12 h[[8]]/(h[[2]] h[[3]]) /;  h[[1]] > 0 ; 
     (* = x1 wL12 ff/(EL1 EL2) *) 
bL2bL1 [ h_] := BL2BL1 /;  h[[1]]==0 ;

Selection of the equations for solving for the parameter values at lattice saturation.

chooseNEandCP::usage = "A routine to determine which sets of equations 
for the neighbor-effect parameters and conditional probabilities should 
be used when r1 = r1Sat." ; 
chooseNEandCP := { 
     If [ type==1, 
          If [ wL1 >0 && wL2 >0 && wL12 >0 && !(wL1==wL2==wL12), 
               NE := NE1 ; NEchoice = "NE1" ] ; 
          If [ wL1==0 && wL2 >0 && wL12 >0, NE := NE2 ; NEchoice = "NE2" ] ; 
          If [ wL1 >0 && wL2==0 && wL12 >0, NE := NE3 ; NEchoice = "NE3" ] ; 
          If [ wL1 >0 && wL2 >0 && wL12==0, NE := NE4 ; NEchoice = "NE4" ] ; 
          If [ wL1==wL2==0 && wL12 >0, NE := NE5 ; NEchoice = "NE5" ] ; 
          If [ wL1==wL2==wL12 || wL1==wL12==0 || wL2==wL12==0, 
               NE := NE6 ; NEchoice = "NE6" ] ; 
          If [ !(wL1==wL12==0 || wL2==wL12==0), 
               CP := CP1 ; CPchoice = "CP1" ] ; 
          If [ wL1==wL12==0 && wL2 >0, CP := CP2 ; CPchoice = "CP2" ] ; 
          If [ wL2==wL12==0 && wL1 >0, CP := CP3 ; CPchoice = "CP3" ] ; 
          If [ wL1==wL2==wL12==0, CP := CP4 ; CPchoice = "CP4" ] ] ; 
     If [ type==2 || type==3, NE := NE7 ; NEchoice = "NE7" ; 
          CP := CP5 ; CPchoice = "CP5" ] ; 
     If [ NEchoice=="NE6" || NEchoice=="NE7", 
          lastPoint = True, lastPoint = False ] ; 
     r1 = r1Sat ; 
     If [ flag, 
          If [ lastPoint, NE ; CP ; finalr2L2fChoose, 
               {EL1,EL2,finalr2L2f} = {1,1,0} ] , 
          NE ; CP ; finalr2L2fChoose ; 
          If [ lastPoint==False, checkLimits ] ]  ; 
     PrependTo [ abscissa, 0 ] ; 
     AppendTo [ abscissa, r1Sat ]  } ;

finalr2L2fChoose::usage = "A routine to determine the value of
r2/L2f when r1 = r1Sat." ; 
finalr2L2fChoose := { 
     If [ type==1, 
          If [ wL1==0 && !wL12==0 && nL2==1 && r2 < r1,
               r2b = r2 (nL1+1) ; finalL2f = r2b/(kL2 wL12^2 (1-r2b)) ; 
               finalr2L2f = r2/finalL2f, 
               finalr2L2f = 0 ] ] ; 
     If [ type==2, finalr2L2f = r2/FreeL2 ] ; 
     If [ type==3, finalr2L2f = r2/(totalL2-r2 polymer) ] } ;

Expressions for the neighbor-effect parameters when r1 is at its maximum value

NE1::usage = "NE1 gives equations for the neighbor-effect parameters
at saturation (used for type 1 isotherms when wL1, wL2 and wL12 are all
nonzero but not all the same)." ; 
NE1 := { 
      {EL1,EL2} = {eL1,eL2}/.FindRoot [ { 
            eL1 == (wL1 r1 eL2 + wL12 r2 eL1)/(r1 eL2 + r2 eL1) , 
            eL2 == (wL2 r2 eL1 + wL12 r1 eL2)/(r1 eL2 + r2 eL1)  } , 
            {eL1,EL1}, {eL2,EL2}, AccuracyGoal->7, 
            MaxIterations ->200, WorkingPrecision ->18 ]  } ;

NE2::usage = "NE2 gives equations for the neighbor-effect parameters
at saturation for type 1 isotherms when wL1 = 0 and wL2 and wL12
are > zero)." ; 
NE2 := {
     If [ r1==r2 || r1 > r2, EL1 = 0 ; EL2 = wL12, 
           {EL1,EL2} = {eL1,eL2}/.FindRoot [ { 
                eL1 == (wL12 r2 eL1)/(r1 eL2 + r2 eL1) , 
                eL2 == (wL2 r2 eL1 + wL12 r1 eL2)/(r1 eL2 + r2 eL1)  } , 
                {eL1,EL1}, {eL2,EL2}, AccuracyGoal->7, 
                MaxIterations ->200, WorkingPrecision ->18 ] ]  } ;

NE3::usage = "NE3 gives equations for the neighbor-effect parameters
for type 1 isotherms at saturation when wL2 = 0 and wL1 and wL12
are > zero)." ; 
NE3 := {
     If [ r1==r2 || r1 < r2, EL1 = wL12 ; EL2 = 0, 
           {EL1,EL2} = {eL1,eL2}/.FindRoot [ { 
                eL1 == (wL1 r1 eL2 + wL12 r2 eL1)/(r1 eL2 + r2 eL1) , 
                eL2 == (wL12 r1 eL2)/(r1 eL2 + r2 eL1)  } , 
                {eL1,EL1}, {eL2,EL2}, AccuracyGoal->7, 
                MaxIterations ->200, WorkingPrecision ->18 ] ]  } ;

NE4::usage = "NE4 gives equations for the neighbor-effect parameters
for type 1 isotherms at saturation when wL12 = 0 and wL1 and wL2
are > zero)." ; 
NE4 := { 
     {EL1,EL2} = {eL1,eL2}/.FindRoot [ { 
          eL1 == (wL1 r1 eL2)/(r1 eL2 + r2 eL1) , 
          eL2 == (wL2 r2 eL1)/(r1 eL2 + r2 eL1)  } , 
          {eL1,EL1}, {eL2,EL2}, AccuracyGoal ->7, 
          MaxIterations ->200, WorkingPrecision ->18 ]  } ;

NE5::usage = "NE5 gives equations for the neighbor-effect parameters
for type 1 isotherms at saturation (used when wL1 = wL2 = 0 and
wL12 > zero)." ; 
NE5 := {
     If [ r1==r2, 
          {EL1,EL2} = {eL1,eL2}/.FindRoot [ { 
               eL1 == (wL12 r2 eL1)/(r1 eL2 + r2 eL1) , 
               eL2 == (wL12 r1 eL2)/(r1 eL2 + r2 eL1)  } , 
               {eL1,EL1}, {eL2,EL2}, AccuracyGoal ->7, 
               MaxIterations ->200, WorkingPrecision ->18 ] , 
          If [ r1 > r2, EL1 = 0 ; EL2 = wL12, EL1 = wL12 ; EL2 = 0 ] ]  } ;

NE6::usage = "NE6 gives equations for the neighbor-effect parameters
for type 1 isotherms at saturation (used when wL1 = wL2 = wL12, or
when wL12 = 0 and either wL1 or wL2 is zero)." ; 
NE6 := EL1 = EL2 = wL12 ;

NE7::usage = "NE7 gives equations for the neighbor-effect parameters
and r2 for type 2 or 3 isotherms at saturation." ; 
NE7 := { EL1 = wL1 ; EL2 = wL12 ; 
        If [ wL1==0 && !wL12==0 && nL2==1, 
             t = totalL2 ; p = polymer/(nL1+1) ; kd = 1/(kL2 wL12^2) ; 
             If [ type==2, r2b = freeL2/(kd + freeL2) ] ;  
             If [ type==3, 
                  r2b = (t+p+kd - Sqrt [ (t+p+kd)^2 - 4 p t ] )/(2 p)  ] ; 
             r2 = r2b/(nL1+1), 
             r2 = 0 ]  } ;

Expressions for conditional probabilities when r1 is at its maximum value

fMin::usage = "fMin is the value of f when the lattice is saturated." ; 
BL1F::usage = "BL1F is the value of bL1f when the lattice is saturated." ; 
BL2F::usage = "BL2F is the value of bL2f when the lattice is saturated." ; 
FBL1::usage = "FBL1 is the value of fbL1 when the lattice is saturated." ; 
FBL2::usage = "FBL2 is the value of fbL2 when the lattice is saturated." ; 
BL1BL1::usage = "BL1BL1 is the value of bL1bL1 when the lattice is
saturated." ; 
BL2BL2::usage = "BL2BL2 is the value of bL2bL2 when the lattice is
saturated." ; 
BL1BL2::usage = "BL1BL2 is the value of bL1bL2 when the lattice is
saturated." ; 
BL2BL1::usage = "BL2BL1 is the value of bL2bL1 when the lattice is
saturated." ;

CP1::usage = "CP1 gives expressions for the conditional probabilities
for type 1 isotherms at saturation when no more than two of wL1, wL2,
and wL12 are zero (and if two are zero, wL12 is > zero)." ; 
CP1 := { BL1F := fMin EL2 / (fMin EL1 EL2 + r1 EL2 + r2 EL1) ; 
     BL2F := fMin EL1 / (fMin EL1 EL2 + r1 EL2 + r2 EL1) ; 
     FBL1 := (r1 EL2) / (fMin EL1 EL2+ r1 EL2 + r2 EL1) ; 
     FBL2 := (r2 EL1) / (fMin EL1 EL2+ r1 EL2 + r2 EL1) ; 
     If [ wL1==0, BL1BL1 = 0, 
          If [ wL2==0 && wL12 >0 && (r2 >r1 || r2==r1), BL1BL1 = 0, 
               BL1BL1 := 1 - BL1F - BL1BL2 ] ] ; 
     If [ wL2==0, BL2BL2 = 0, 
          If [ wL1==0 && wL12 >0 && (r1 >r2 || r1==r2), BL2BL2 = 0, 
               BL2BL2 := 1 - BL2F - BL2BL1 ] ] ; 
     BL1BL2 := r2 wL12 / (fMin EL1 EL2 + r1 EL2 + r2 EL1) ; 
     BL2BL1 := r1 wL12 / (fMin EL1 EL2 + r1 EL2 + r2 EL1)  } ;

CP2::usage = "CP2 gives the values of the conditional probabilities
for type 1 isotherms at saturation when wL1 = wL12 = 0 and wL2 > 0." ; 
CP2 := { BL1F = 1 ; BL2F = 0 ; FBL1 = 1 ; FBL2 = 0 ; 
            BL1BL1 = 0 ; BL2BL2 = 1 ; BL1BL2 = 0 ; BL2BL1 = 0  } ;

CP3::usage = "CP3 gives the values of the conditional probabilities
for type 1 isotherms at saturation when wL2 = wL12 = 0 and wL1 > 0." ; 
CP3 := { BL1F = 0 ; BL2F = 1 ; FBL1 = 0 ; FBL2 = 1 ; 
            BL1BL1 = 1 ; BL2BL2 = 0 ; BL1BL2 = 0 ; BL2BL1 = 0  } ;

CP4::usage = "CP4 gives the values of conditional probabilities
for type 1 isotherms at saturation when wL1, wL2, and wL12 are all zero." ; 
CP4 := { BL1F = 1 ; BL2F = 1 ; FBL1 = r1/(r1+r2) ; FBL2 = r2/(r1+r2) ; 
            BL1BL1 = 0 ; BL2BL2 = 0 ; BL1BL2 = 0 ; BL2BL1 = 0  } ;

CP5::usage = "CP5 gives the values of the conditional probabilities
for type 2 or 3 isotherms at saturation." ; 
CP5 := { FBL1 = 1 ; FBL2 = 0 ; BL2BL2 = 0 ; 
            If [ wL12==0, BL2BL1 = 0 ; BL2F = 1, 
                 BL2BL1 = 1 ; BL2F = 0 ] ; 
            If [ wL1==0, BL1BL1 = 0 ; 
                 If [ !wL12==0 && nL2==1, 
                      BL1F = 1-r2b ; BL1BL2 = r2b,
                      BL1F = 1 ; BL1BL2 = 0 ], 
                 BL1BL1 = 1 ; BL1F = 0 ; BL1BL2 = 0 ] } ;

Define plotting routines

makePlot::usage = "makePlot is a neighbor-effect parameter plotting
routine." ; 
makePlot [ j_, yOrigin_, yLabel_, graphLabel_ ] := { 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, 
          If [ j==1, 
               If [ r1ZeroError, printWarning2, printWarning1 ], 
               printWarning2 ] ; 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,j]]}, 
               {i,points-dropPoints} ] ; 
          If [ lastPoint || j==1 || j==8, 
               AppendTo [ plotTable, {abscissa[[points+1]], 
                    ordinate[[points+1,j]]} ] ] , 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,j]]}, {i,points+1} ] ] ; 
     If [ r1ZeroError==False, 
          curvePlot = ListPlot [ plotTable, PlotJoined->True, 
               AxesLabel->{" r1",yLabel}, PlotLabel->graphLabel, 
               DisplayFunction->Identity, 
               If [ j==8, 
                    PlotRange->{{-0.025 r1Sat,1.025 r1Sat},{-0.025,1.025}}, 
                    PlotRange->All ] ] ;
          pointPlot = ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
               DisplayFunction->Identity, 
               If [ j==8, 
                    PlotRange->{{-0.025 r1Sat,1.025 r1Sat},{-0.025,1.025}}, 
                    PlotRange->All ] ] ;
          Show [ curvePlot, pointPlot, DisplayFunction->$DisplayFunction, 
               If [ versionIsLessThan2, Axes->{0,yOrigin}, 
                    AxesOrigin->{0,yOrigin} ] ] ] ; 
     N [ plotTable ] }

makeCPplot::usage = "makeCPplot is a plotting routine for conditional
probabilities." ; 
makeCPplot [ cpFunction_, label_ ] := { 
     cpTable = Map [ cpFunction , ordinate ] ; 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, printWarning2 ; 
          plotTable = Table [ {abscissa[[j]], cpTable[[j]]}, 
               {j,points-dropPoints} ] ; 
          If [ lastPoint, 
               AppendTo [ plotTable, 
                    {abscissa[[points+1]], cpTable[[points+1]]} ] ], 
          plotTable = Table [ {abscissa[[j]], cpTable[[j]]}, { j, points+1} ] ] ; 
     If [ r1ZeroError==False, 
          curvePlot := ListPlot [ plotTable, PlotJoined->True, 
               AxesLabel->{" r1", label}, 
               DisplayFunction->Identity, 
               PlotRange->{{-0.025 r1Sat,1.025 r1Sat},{-0.025,1.025}} ] ; 
          pointPlot := ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
               DisplayFunction->Identity, 
               PlotRange->{{-0.025 r1Sat,1.025 r1Sat},{-0.025,1.025}} ] ; 
          Show [ curvePlot, pointPlot, DisplayFunction->$DisplayFunction] ] ; 
     N [ plotTable ] }