The Binding of Symmetric Ligands of Two
Kinds to an Infinite Isotropic Linear Lattice
by Alan R. Wolfe
Department of Biopharmaceutical Sciences
University of California
San Francisco, CA  94143-0446

This is the initialization notebook for the case of competitive binding of 2 kinds of symmetic ligand to an isotropic lattice.  This version DOES NOT allow comparison of experimental and theoretical data.  The notebook needs to be evaluated only once per Mathematica  session.  After opening this notebook, execute the "Kernel-> Evaluation-> Evaluate Initialization" command.  Then close this notebook (without quitting Mathematica ) and open the corresponding 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 ] ;

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)]; latticeCoverage := r1 nL1 + r2 nL2;

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 == 3, 
          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/2) (1 - r2 (nL2 - 1)), 
              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 == 1 || type == 2, 
          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==3, If [ wL1==wL2==wL12, 
          ne := ne2 ; neChoice = "ne2", ne := ne1 ; neChoice = "ne1" ] ] ; 
     If [ type==1 || type==2, 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 'totalL2' (for type 1" ] ; 
Print [ "isotherms) or 'freeL2' (for type 2 isotherms)." ] } ;

ne1::usage = "ne1 gives equations for the neighbor-effect parameters for an unsaturated lattice (used for type 3 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 3 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 3 isotherms when wL1 = wL2 = wL12." ;

ne2 := 
      EL1 = EL2 = 
          (1/2) (
              1 - x1 - x2 + 
                (Sqrt[(((1 - x1 - x2))^2) + 4 (x1 + x2) wL1]));

ne3::usage = "ne3 gives equations for the neighbor-effect parameters for an unsaturated lattice (used for type 1 or 2 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 1 or 2 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 1 or 2 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 1 isotherms)." ;

neL2a := {{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]};

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 2 isotherms)." ;

neL2b := {{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))};

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==3, 
          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==1 || type==2, 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, finalr2L2f = (r2/(totalL2 - r2 polymer))]; 
        If[type == 2, finalr2L2f = (r2/FreeL2)]; 
        If[type == 3, 
          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]]};

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 3 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 3 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 3 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 3 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 3 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 3 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 1 or 2 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 == 1, 
            r2b = ((t + p + kd - (Sqrt[(((t + p + kd))^2) - 4 p t]))/(2 p
                  ))]; If[type == 2, r2b = (freeL2/(kd + freeL2))]; 
          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 3 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 3 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 3 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 3 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 1 or 2 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,points2-dropPoints} ] ; 
          If [ lastPoint || j==1 || j==8, 
               AppendTo [ plotTable, {abscissa[[points2+1]], 
                    ordinate[[points2+1,j]]} ] ] , 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,j]]}, {i,points2+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, 
               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,points2-dropPoints} ] ; 
          If [ lastPoint, 
               AppendTo [ plotTable, 
                    {abscissa[[points2+1]], cpTable[[points2+1]]} ] ], 
          plotTable = Table [ {abscissa[[j]], cpTable[[j]]}, { j, points2+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 ] }