The Binding of Asymmetric Ligands to
an Infinite Anisotropic 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 anisotropic lattice, asymmetric ligand case with two binding modes (n1 need not equal n2).  It 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.  It is best to quit Mathematica  and start a new session before evaluating a different initialization notebook.

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 values to be calculated for each of the parameters displayed in the output.  Since the values at r = 0 are not calculated, it equals the number of intervals on the x-axis of the plots." ; 
k1::usage = "k1 is the binding constant for interactions between the lattice and ligands with their heads to the left." ; 
k2::usage = "k2 is the binding constant for interactions between the lattice and ligands with their heads to the right." ; 
n1::usage = "n1 is the number of lattice residues occupied per b1L-oriented bound ligand." ; 
n2::usage = "n2 is the number of lattice residues occupied per b2L-oriented bound ligand." ; 
n::usage = "n is the effective binding site size for the reference symmetric ligand binding Scatchard plot (n = n1 = n2 except when both w12 and w21 and one of w11 or w22 are zero, in which case n = n1 + 1/2)." ; 
w::usage = "w is the value of the cooperativity parameter for the reference symmetric ligand Scatchard plot." ; 
w11::usage = "w11 is the cooperativity parameter for head-head (b1R-b1L) interactions between asymmetric ligands bound to an anisotropic lattice." ; 
w22::usage = "w22 is the cooperativity parameter for tail-tail (b2R-b2L) interactions between asymmetric ligands bound to an anisotropic lattice." ; 
w12::usage = "w12 is the cooperativity parameter for head-tail (b1R-b2L) interactions between asymmetric ligands bound to an anisotropic lattice." ; 
w21::usage = "w21 is the cooperativity parameter for tail-head (b2R-b1L) interactions between asymmetric ligands bound to an anisotropic lattice." ; 
r::usage = "r is the ratio of bound ligands to lattice residues (plotted on the abscissa of a Scatchard plot)." ; 
f::usage = "f is the fraction of lattice residues that are free (used when solving equations)." ; 
F::usage = "F is the fraction of lattice residues that are free (used when not solving equations)." ; 
x::usage = "x is r/f." ; 
X::usage = "X is r/F." ; 
Q::usage = "Q is the fraction of all bound ligands with their heads to the left (b1L-oriented)." ; 
rLf::usage = "rLf is r divided by the free ligand concentration (plotted on the ordinate of a Scatchard plot)." ; 
E1L::usage = "E1L is the neighbor-effect parameter for the left end of a binding site on an anisotropic lattice to be bound by a b1L-oriented ligand." ; 
E1R::usage = "E1R is the neighbor-effect parameter for the right end of a binding site on an anisotropic lattice to be bound by a b2L-oriented ligand." ; 
E2L::usage = "E2L is the neighbor-effect parameter for the left end of a binding site on an anisotropic lattice to be bound by a b2L-oriented ligand." ; 
E2R::usage = "E2R is the neighbor-effect parameter for the right end of a binding site on an anisotropic lattice to be bound by a b1L-oriented ligand." ; 
rSat::usage = "rSat is the limiting value of r, i.e., the ratio of ligands to lattice residues when the lattice is saturated." ; 
deltaA::usage = "deltaA is the quantity
|2(k2 w12 - k1 w21)-(k2 - k1)(w11 + w22)|/(k1 + k2)  -  |w11 - w22|.
(This is a dimensionless version of the parameter that was defined in the JMB paper.)  When deltaA is positive, zero, or negative, the asymmetric ligand Scatchard plot will be above, the same as, or below its reference symmetric ligand Scatchard plot, respectively.  (The reference plot is the unique symmetric ligand Scatchard plot with the same intercepts and initial slope.)  On an isotropic lattice, deltaA is always less than or equal to zero." ;

\(x := r\/f; 

      rLfRef := 
        1\/2\ \((k1 + k2)\)\ \((1 - n\ r)\)\ e\^2\ 
          \((1 + xRef\/e)\)\^\(1 - n\); \)

Def1::usage = "Definitions for the case n1 = n2." ;

Def1 := {f := 1 - n1 r; F := 1 - n1 r; 
        rLf := Re[
            (1/2) F (ff^(n1 - 1)) (k1 E1L E2R + k2 E1R E2L)]; 
        deltaA1 := 
          Abs[2 (k2 w12 - k1 w21) - (k2 - k1) (w11 + w22)]
                /(k1 + k2) - Abs[w11 - w22]; 
        If[k1 > 0 && k2 > 0 && Abs[deltaA1] ≥ (1/(10^12)), 
          deltaA := deltaA1, deltaA = 0]};

Def2::usage = "Definitions for the case n1 not equal to n2." ;

Def2 := {f := 1 - n1 q r - n2 (1 - q) r; 
        F := 1 - n1 Q r - n2 (1 - Q) r; 
        rLf := Re[
            (1/2) F (
                (ff^(n1 - 1)) k1 E1L E2R + (ff^(n2 - 1)) k2 E1R E2L)
                ]; If[k1 > 0 && k2 > 0, deltaA = "undefined", 
          deltaA = 0]};

Calculation of the initial slope

initialSlope := 
      ((k1 k2 (w11 + w22 - w12 - w21))/(k1 + k2)) + 
        k1 (w21 - n1 - 0.5) + k2 (w12 - n2 - 0.5);

Determination of the maximum value of r

rSatDetermine::usage = "A routine to determine the ratio of bound ligands to lattice residues when the lattice has reached saturation." ;

rSatDetermine := {nmin = Min[n1, n2]; rSat = (1/nmin); 
        If[n1 < n2 && w21 == 0, rSat = (1/(n1 + 1))]; 
        If[n2 < n1 && w12 == 0, rSat = (1/(n2 + 1))]; 
        If[w11 == w22 == w12 == w21 == 0 && k1 > 0 && k2 > 0, 
          rSat = (1/(nmin + 1))]; 
        If[k1 == 0, If[w12 == 0, rSat = (1/(n2 + 1)), rSat = (1/n2)]]; 
        If[k2 == 0, If[w21 == 0, rSat = (1/(n1 + 1)), rSat = (1/n1)]]; 
        If[w12 == w21 == 0 && Xor[w11 == 0, w22 == 0] && k1 > 0 && k2 > 0, 
          rSat = Max[(2/(n1 + n2 + 1)), (1/(n1 + 1)), (1/(n2 + 1))]]; 
        If[w12 == w21 == 0 && w11 > 0 && w22 > 0 && k1 > 0 && k2 > 0, 
          rSat = Max[(2/(n1 + n2)), (1/(n1 + 1)), (1/(n2 + 1))]]; 
        If[w12 == 0 && w21 > 0 && w11 > 0 && w22 > 0 && k1 > 0 && k2 > 0, 
          rSat = Max[(2/(n1 + n2)), (1/n1), (1/(n2 + 1))]]; 
        If[w21 == 0 && w12 > 0 && w11 > 0 && w22 > 0 && k1 > 0 && k2 > 0, 
          rSat = Max[(2/(n1 + n2)), (1/n2), (1/(n1 + 1))]]};

Equations for calculating values of neighbor-effect parameters and q when r is below its maximum value

neChoose::usage = "A routine to determine which set of equations for neighbor-effect parameters should be used when r < rSat, to decide whether the reference symmetric ligand Scatchard plot needs to be calculated, and to determine the parameter values for the reference Scatchard plot." ; 
neChoose := { 
     flag = False ; setLimits ; 
     w = k1 k2 (w11+w22-w12-w21)/(k1+k2)^2 + 
          (k1 w21 + k2 w12)/(k1 + k2) ; 
     If [ Abs [deltaA] <10^(-12) || k1==0 || k2==0 || !n1==n2, 
          needReference = False, 
          needReference = True ; 
          n = n1 ; 
          If [ w12==0 && w21==0 && 
               Xor [ w11==0, w22==0 ], w = w + 0.5 ; n = n + 0.5 ] ] ; 
     If [ n1==n2 && k1 >0 && k2 >0, 
          If [ w11==w22==w12==w21==1, ne := ne1a ; neChoice = "ne1a", 
               If [ w11==w22, 
                    If [ Abs [ deltaA ] >10^(-12), ne := ne2a ; neChoice = "ne2a", 
                         ne := ne2b ; neChoice = "ne2b" ], 
                    If [ k1==k2 && w12==w21, ne := ne1b ; neChoice = "ne1b", 
                         If [ Abs [ deltaA ] >10^(-12), ne := ne1 ; neChoice = "ne1", 
                              ne := ne1c ; neChoice = "ne1c" ] ] ] ] ] ; 
     If [ k1==0, ne := ne3 ; neChoice = "ne3"  ] ; 
     If [ k2==0, ne := ne4 ; neChoice = "ne4"  ] ; 
     If [ !n1==n2 && k1 >0 && k2 >0, 
          If [ w11==w22==w12==w21==1, 
               If [ n1 >n2, ne := ne7a ; neChoice = "ne7a", 
                    ne := ne7b ; neChoice = "ne7b" ], 
               If [ w11==w22, ne := ne6 ; neChoice = "ne6", 
                    ne := ne5 ; neChoice = "ne5" ]  ]  ] ; 
     If [ neChoice=="ne1" || neChoice=="ne1b" || neChoice=="ne2a" || 
          neChoice=="ne5" || neChoice=="ne6" || neChoice=="ne7a" || 
          neChoice=="ne7b", needCheck = True, needCheck = False ]  } ;

setLimits::usage = "setLimits is a routine that determines the largest and smallest physically possible values for the neighbor-effect parameters." ; 
setLimits := { 
     neTol = 0.0000001 ; 
     E1Lmin = Min [ 1, w11, w21 ] - neTol ; 
     E1Lmax = Max [ 1, w11, w21 ] + neTol ; 
     E1Rmin = Min [ 1, w11, w12 ] - neTol ; 
     E1Rmax = Max [ 1, w11, w12 ] + neTol ; 
     E2Lmin = Min [ 1, w22, w12 ] - neTol ; 
     E2Lmax = Max [ 1, w22, w12 ] + neTol ; 
     E2Rmin = Min [ 1, w22, w21 ] - neTol ; 
     E2Rmax = Max [ 1, w22, w21 ] + 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, Q and the 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) || (r < rError), 
          If [ Head [ E1L ] == Complex || Head [ E1R ] == Complex || 
               Head [ E2L ] == Complex || Head [ E2R ] == Complex || 
               Head [ Q ] == Complex, 
               {E1L,E1R,E2L,E2R,Q} = Map [ Re, {E1L,E1R,E2L,E2R,Q} ] ; 
               rError = r ; 
               If [ flag==False, needPrintWarning1 = True ] ; 
               flag = True ]  ] ; 
     If [ (flag==False) || (r < rError), 
          If [ E1L<E1Lmin || E1L>E1Lmax || E1R<E1Rmin || E1R>E1Rmax || 
               E2L<E2Lmin || E2L>E2Lmax || E2R<E2Rmin || E2R>E2Rmax || 
               Q<0 || Q>1 || F<0 || F>1, 
               rError = r ; 
               If [ flag==False, needPrintWarning1 = True ] ; 
               flag = True ]  ] ; 
     If [ ((flag==False) || (r < rError)) && r < rSat, 
          X = r/F ; fF = ff ; b1LfTest = fF/E1L ; b2LfTest = fF/E2L ; 
          b1RfTest = fF/E1R ; b2RfTest = fF/E2R ; 
          fb1Ltest = Q X fF/E1L ; fb1Rtest = (1-Q) X fF/E1R ; 
          fb2Ltest = (1-Q) X fF/E2L ; fb2Rtest = Q X fF/E2R ; 
          b1Lb1Rtest = (1-Q) X w11 fF/(E1L E1R) ; 
          b1Rb1Ltest = Q X w11 fF/(E1R E1L) ; 
          b2Lb2Rtest = Q X w22 fF/(E2L E2R) ; 
          b2Rb2Ltest = (1-Q) X w22 fF/(E2R E2L) ; 
          b1Lb2Rtest = Q X w21 fF/(E1L E2R) ; 
          b1Rb2Ltest = (1-Q) X w12 fF/(E1R E2L) ] ; 
     If [ flag==False && r==rSat, 
          fF = 0 ; b1LfTest = B1LF ; b2LfTest = B2LF ; 
          b1RfTest = B1RF ; b2RfTest = B2RF ; 
          fb1Ltest = FB1L ; fb1Rtest = FB1R ; 
          fb2Ltest = FB2L ; fb2Rtest = FB2R ; 
          b1Lb1Rtest = B1LB1R ; b1Rb1Ltest = B1RB1L ; 
          b2Lb2Rtest = B2LB2R ; b2Rb2Ltest = B2RB2L ; 
          b1Lb2Rtest = B1LB2R ; b1Rb2Ltest = B1RB2L ] ; 
     If [ (flag==False) || (r < rError), 
          If [ fF < cpMin || fF > cpMax || 
               b1LfTest < cpMin || b1LfTest > cpMax1 || 
               b2LfTest < cpMin || b2LfTest > cpMax || 
               b1RfTest < cpMin || b1RfTest > cpMax || 
               b2RfTest < cpMin || b2RfTest > cpMax || 
               fb1Ltest < cpMin || fb1Ltest > cpMax || 
               fb2Ltest < cpMin || fb2Ltest > cpMax || 
               fb1Rtest < cpMin || fb1Rtest > cpMax || 
               fb2Rtest < cpMin || fb2Rtest > cpMax || 
               b1Lb1Rtest < cpMin || b1Lb1Rtest > cpMax || 
               b1Rb1Ltest < cpMin || b1Rb1Ltest > cpMax || 
               b2Lb2Rtest < cpMin || b2Lb2Rtest > cpMax || 
               b2Rb2Ltest < cpMin || b2Rb2Ltest > cpMax || 
               b1Lb2Rtest < cpMin || b1Lb2Rtest > cpMax || 
               b1Rb2Ltest < cpMin || b1Rb2Ltest > cpMax, 
               rError = r ; 
               If [ flag==False, needPrintWarning1 = True ] ; 
               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 r = ", rError//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." ] ; 
Print [ " " ] } ;

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

ne1::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 does not = w22, deltaA does not = 0, and n1 = n2." ;

ne1 := {
        Z = ((k1 E1L E2R + k2 E1R E2L)/x); {E1L, E1R, E2L, 
            E2R} = {e1l, e1r, e2l, e2r} /. 
            FindRoot[{
                e1l (z + k2 e2l + k1 e1l) == 
                  z + w11 k2 e2l + w21 k1 e1l, 
                e1r (z + k1 e2r + k2 e1r) == 
                  z + w11 k1 e2r + w12 k2 e1r, 
                e2l (z + k2 e2l + k1 e1l) == 
                  z + w12 k2 e2l + w22 k1 e1l, 
                e2r (z + k1 e2r + k2 e1r) == 
                  z + w21 k1 e2r + w22 k2 e1r, 
                z x == k1 e1l e2r + k2 e1r e2l}, {e1l, E1L}, {e1r, 
                E1R}, {e2l, E2L}, {e2r, E2R}, {z, Z}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]; 
        Q = ((k1 E1L E2R)/(k1 E1L E2R + k2 E1R E2L))};

ne1::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 does not = w22, deltaA does not = 0, and n1 = n2." ; 
(* Alternate Version *)

ne1 := {
        Z = ((k1 E1L E2R + k2 E1R E2L)/x); {E1L, E1R, E2L, 
            E2R} = {e1l, e1r, e2l, e2r} /. 
            FindRoot[{
                e1l == ((z + w11 k2 e2l + w21 k1 e1l)/(z + k2 e2l + 
                      k1 e1l)), 
                e1r (z + k1 e2r + k2 e1r) == 
                  z + w11 k1 e2r + w12 k2 e1r, 
                z (1 - e2l) == 
                  k1 e1l (e2l - w22) + k2 e2l (e2l - w12), 
                z (1 - e2r) + k1 e2r (w21 - e2r) + 
                    k2 e1r (w22 - e2r) == 0, 
                z x == k1 e1l e2r + k2 e1r e2l}, {e1l, E1L}, {e1r, 
                E1R}, {e2l, E2L}, {e2r, E2R}, {z, Z}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]; 
        Q = ((k1 E1L E2R)/(k1 E1L E2R + k2 E1R E2L))};

ne1::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 does not = w22, deltaA does not = 0, and n1 = n2." ; 
(* Alternate Version *)

ne1 := {{E1L, E1R, E2L, E2R} = {e1l, e1r, e2l, e2r} /. 
            FindRoot[{
                (k1 e1l e2r + k2 e1r e2l) (1 - e1l) == 
                  k1 x e1l (e1l - w21) + k2 x e2l (e1l - w11), 
                (k1 e1l e2r + k2 e1r e2l) (1 - e1r) == 
                  k1 x e2r (e1r - w11) + k2 x e1r (e1r - w12), 
                k1 e1l (e2r + x) + k2 e2l (e1r + x) == 
                  ((k1 e1l (e2r + x w22))/e2l) + 
                    k2 (e1r + x w12), 
                k1 e2r (e1l + x) + k2 e1r (e2l + x) == 
                  k1 (e1l + x w21) + 
                    ((k2 e1r (e2l + x w22))/e2r)}, {e1l, E1L}, {e1r, 
                E1R}, {e2l, E2L}, {e2r, E2R}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]; 
        Q = ((k1 E1L E2R)/(k1 E1L E2R + k2 E1R E2L))};

ne1::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 does not = w22, deltaA does not = 0, and n1 = n2." ; 
(* Alternate Version *)

ne1 := {{E1L, E1R, E2L, E2R, Q} = {e1l, e1r, e2l, e2r, q} /. 
            FindRoot[{
                e1l == (((1/x) + ((w11 (1 - q))/e1r) + 
                      ((w21 q)/e2r))/((1/x) + ((1 - q)/e1r) + (q/e2r))), 
                e1r == (((1/x) + ((w11 q)/e1l) + 
                      ((w12 (1 - q))/e2l))/((1/x) + (q/e1l) + 
                      ((1 - q)/e2l))), 
                e2l == (((1/x) + ((w12 (1 - q))/e1r) + 
                      ((w22 q)/e2r))/((1/x) + ((1 - q)/e1r) + (q/e2r))), 
                e2r == (((1/x) + ((w21 q)/e1l) + 
                      ((w22 (1 - q))/e2l))/((1/x) + (q/e1l) + 
                      ((1 - q)/e2l))), 
                k1 (e1l - e2r) == k2 (e1r - e2l)}, {e1l, E1L}, {
                e1r, E1R}, {e2l, E2L}, {e2r, E2R}, {q, Q}, 
              AccuracyGoal -> 7, MaxIterations -> 200, 
              WorkingPrecision -> 17]};

ne1a::usage = "Neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 = w22 = w12 = w21 = 1, and n1 = n2." ; 
ne1a := {E1L,E1R,E2L,E2R,Q} = {1,1,1,1, k1/(k1 + k2) } ;

ne1b::usage = "Equations for neighbor-effect parameters when r < rSat, k1 = k2, n1 = n2, w12 = w21 (the lattice is isotropic), and w11 does not = w22." ;

ne1b := {{E1L, E1R, E2L, E2R, Q} = {e1, e1, e2, e2, 0.5} /. 
            FindRoot[{
                e1 ((2/x) + (1/e1) + (1/e2)) == (2/x) + (w11/e1) + (w12/e2), 
                e2 ((2/x) + (1/e1) + (1/e2)) == (2/x) + (w12/e1) + (w22/e2)}, {
                e1, E1L}, {e2, E2L}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]};

ne1b::usage = "Equations for neighbor-effect parameters when r < rSat, k1 = k2, n1 = n2, w12 = w21 (the lattice is isotropic), and w11 does not = w22." ; 
(* Alternate Version *)

ne1b := {{E1L, E1R, E2L, E2R, Q} = {e1, e1, e2, e2, 0.5} /. 
            FindRoot[{
                e1 == (((2/x) + (w11/e1) + (w12/e2))/((2/x) + (1/e1) + (1/e2))), 
                e2 == (((2/x) + (w12/e1) + (w22/e2))/((2/x) + (1/e1) + 
                      (1/e2)))}, {e1, E1L}, {e2, E2L}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]};

ne1c::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 does not = w22, deltaA = 0, and n1 = n2." ;

ne1c := {e = (1/2) (1 - x + (Sqrt[(((1 - x))^2) + 4 x w])); 
        E1L = (((k1 + k2) 
              ((e^2) + w11 e x))/((k1 + k2) ((e^2) + e x) + 
              k1 x (w11 - w21))); 
        E2L = (((k1 + k2) e - k1 E1L)/k2); 
        E1R = (((k1 + k2) 
              ((e^2) + w11 e x))/((k1 + k2) ((e^2) + e x) + 
              k2 x (w11 - w12))); 
        E2R = (((k1 + k2) e - k2 E1R)/k1); 
        Q = ((k1 E1L E2R)/(k1 E1L E2R + k2 E1R E2L))};

ne2a::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11, w22, w12 and w21 are not all one, w11 = w22, deltaA does not = 0, and n1 = n2." ;

ne2a := {
        Z = ((k1 (E1L^2) + k2 (E1R^2))/x); {E1L, E1R, E2L, 
            E2R} = {e1l, e1r, e1r, e1l} /. 
            FindRoot[{
                e1l (z + k2 e1r + k1 e1l) == 
                  z + w11 k2 e1r + w21 k1 e1l, 
                e1r (z + k1 e1l + k2 e1r) == 
                  z + w11 k1 e1l + w12 k2 e1r, 
                z x == k1 (e1l^2) + k2 (e1r^2)}, {e1l, E1L}, {e1r, E1R}, {z, 
                Z}, AccuracyGoal -> 7, MaxIterations -> 200, 
              WorkingPrecision -> 17]; 
        Q = ((k1 (E1L^2))/(k1 (E1L^2) + k2 (E1R^2)))};

ne2a::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11, w22, w12 and w21 are not all one, w11 = w22, deltaA does not = 0, and n1 = n2." ; 
(* Alternate Version *)

ne2a := {{E1L, E1R, E2L, E2R, Q} = {e1l, e1r, e1r, e1l, q} /. 
            FindRoot[{
                e1l == (((1/x) + ((w11 (1 - q))/e1r) + 
                      ((w21 q)/e1l))/((1/x) + ((1 - q)/e1r) + (q/e1l))), 
                e1r == (((1/x) + ((w11 q)/e1l) + 
                      ((w12 (1 - q))/e1r))/((1/x) + (q/e1l) + 
                      ((1 - q)/e1r))), 
                ((k1 (1 - q) e1l)/e1r) == ((k2 q e1r)/e1l)}, {e1l, 
                E1L}, {e1r, E1R}, {q, Q}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]};

ne2b::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 = w22, w11, w22, w12 and w21 are not all one, deltaA = 0, and n1 = n2." ;

ne2b := {Q = (k1/(k1 + k2)); 
        E1L = E1R = 
            E2L = E2R = 
                (1/2) (
                    1 - X + 
                      ((((((1 - X))^2) + 
                            ((4 X (k2 w11 + k1 w21))/(k1 + k2)))
                          )^0.5))};

ne3::usage = "Equation for neighbor-effect parameters when r < rSat and k1 = zero." ;

ne3 := {
        E1R = E2L = 
            (1/2) (1 - X + ((((((1 - X))^2) + 4 X w12))^0.5)); 
        Q = 0};

ne4::usage = "Equation for neighbor-effect parameters when r < rSat and k2 = zero." ;

ne4 := {
        E1L = E2R = 
            (1/2) (1 - X + ((((((1 - X))^2) + 4 X w21))^0.5)); 
        Q = 1};

ne5::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 does not = w22, and n1 does not = n2." ;

ne5 := {{E1L, E1R, E2L, E2R, Q} = {e1l, e1r, e2l, e2r, q} /. 
            FindRoot[{
                e1l (z + ((1 - q)/e1r) + (q/e2r)) == 
                  z + ((w11 (1 - q))/e1r) + ((w21 q)/e2r), 
                e1r (z + (q/e1l) + ((1 - q)/e2l)) == 
                  z + ((w11 q)/e1l) + ((w12 (1 - q))/e2l), 
                e2l (z + ((1 - q)/e1r) + (q/e2r)) == 
                  z + ((w12 (1 - q))/e1r) + ((w22 q)/e2r), 
                e2r (z + (q/e1l) + ((1 - q)/e2l)) == 
                  z + ((w21 q)/e1l) + ((w22 (1 - q))/e2l), 
                k1 (e1l - e2r) 
                    (((1 + (q/(z e1l)) + ((1 - q)/(z e2l))))^(n2 - n1
                        )) == k2 (e1r - e2l), 
                z == (1/r) - n1 q - n2 (1 - q)}, {e1l, E1L}, {e1r, 
                E1R}, {e2l, E2L}, {e2r, E2R}, {q, Q}, {z, (1/X)}, 
              AccuracyGoal -> 7, MaxIterations -> 200, 
              WorkingPrecision -> 17]};

ne5::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 does not = w22, and n1 does not = n2." ; (* Alternate Version *)

ne5 := {{E1L, E1R, E2L, E2R, Q} = {e1l, e1r, e2l, e2r, q} /. 
            FindRoot[{
                e1l == (((1/x) + ((w11 (1 - q))/e1r) + 
                      ((w21 q)/e2r))/((1/x) + ((1 - q)/e1r) + (q/e2r))), 
                e1r == (((1/x) + ((w11 q)/e1l) + 
                      ((w12 (1 - q))/e2l))/((1/x) + (q/e1l) + 
                      ((1 - q)/e2l))), 
                e2l == (((1/x) + ((w12 (1 - q))/e1r) + 
                      ((w22 q)/e2r))/((1/x) + ((1 - q)/e1r) + (q/e2r))), 
                e2r == (((1/x) + ((w21 q)/e1l) + 
                      ((w22 (1 - q))/e2l))/((1/x) + (q/e1l) + 
                      ((1 - q)/e2l))), 
                k1 (e1l - e2r) 
                    (((1 + ((q x)/e1l) + (((1 - q) x)/e2l)))^(n2 - 
                        n1)) == k2 (e1r - e2l)}, {e1l, E1L}, {e1r, 
                E1R}, {e2l, E2L}, {e2r, E2R}, {q, Q}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]};

ne6::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 = w22, w11, w22, w12 and w21 are not all one, and n1 does not = n2." ;

ne6 := {{E1L, E1R, E2L, E2R, Q} = {e1l, e1r, e1r, e1l, q} /. 
            FindRoot[{
                e1l (z + ((1 - q)/e1r) + (q/e1l)) == 
                  z + ((w11 (1 - q))/e1r) + ((w21 q)/e1l), 
                e1r (z + (q/e1l) + ((1 - q)/e1r)) == 
                  z + ((w11 q)/e1l) + ((w12 (1 - q))/e1r), 
                ((k1 (1 - q) e1l 
                      (((1 + (q/(z e1l)) + ((1 - q)/(z e1r))))^(n2 - 
                          n1)))/e1r) == ((k2 q e1r)/e1l), 
                z == (1/r) - n1 q - n2 (1 - q)}, {e1l, E1L}, {e1r, 
                E1R}, {q, Q}, {z, (1/X)}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]};

ne6::usage = "Equations for neighbor-effect parameters when r < rSat, both k1 and k2 are > zero, w11 = w22, w11, w22, w12 and w21 are not all one, and n1 does not = n2." ; 
(* Alternate Version *)

ne6 := {{E1L, E1R, E2L, E2R, Q} = {e1l, e1r, e1r, e1l, q} /. 
            FindRoot[{
                e1l == (((1/x) + ((w11 (1 - q))/e1r) + 
                      ((w21 q)/e1l))/((1/x) + ((1 - q)/e1r) + (q/e1l))), 
                e1r == (((1/x) + ((w11 q)/e1l) + 
                      ((w12 (1 - q))/e1r))/((1/x) + (q/e1l) + 
                      ((1 - q)/e1r))), 
                ((k1 (1 - q) e1l 
                      (((1 + ((q x)/e1l) + (((1 - q) x)/e1r))
                          )^(n2 - n1)))/e1r) == ((k2 q e1r)/e1l)}, {
                e1l, E1L}, {e1r, E1R}, {q, Q}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]};

ne7a::usage = "Equations for neighbor-effect parameters when r < max, both k1 and k2 are > zero, w11 = w22 = w12 = w21 = 1, and n1 > n2.  The equivalent equations 7a and 7b are needed because of peculiarities of the Mathematica  function FindRoot." ;

ne7a := {{E1L, E1R, E2L, E2R, Q} = {1, 1, 1, 1, q} /. 
            FindRoot[
              k1 (1 - q) (((1 + x))^(n2 - n1)) == k2 q, {q, Q}, 
              AccuracyGoal -> 12, MaxIterations -> 200, 
              WorkingPrecision -> 17]};

ne7b::usage = "Equations for neighbor-effect parameters when r < max, both k1 and k2 are > zero, w11 = w22 = w12 = w21 = 1, and n1 < n2." ;

ne7b := {{E1L, E1R, E2L, E2R, Q} = {1, 1, 1, 1, q} /. 
            FindRoot[
              k1 (1 - q) == k2 q (((1 + x))^(n1 - n2)), {q, Q}, 
              AccuracyGoal -> 12, MaxIterations -> 200, 
              WorkingPrecision -> 17]};

startPrint::usage = "startPrint starts the printout of the Scatchard plot." ; 
startPrint := { 
     Print [ "r           r/Lf" ] ; 
     Print [ "0          ", (k1+k2)/2//N ] } ;

makeReference::usage = "makeReference defines and maps a function onto the abscissa array to create the reference symmetric ligand Scatchard plot (it checks first to see if this is needed).  It also adds some initial and terminal values to arrays and prints the last line." ; 
makeReference := { 
     PrependTo [ abscissa, 0 ] ; 
     If [ needReference, 
          ordinate2 = Map [ Function [ xx, r = xx ; 
               xRef = r/(1 - n r) ; 
               e = (1-xRef + Sqrt [(1-xRef)^2 + 4 xRef w] )/2 ; 
               rLfRef ], abscissa ] ; 
          AppendTo [ ordinate2, 0 ]  ] ; 
     AppendTo [ abscissa, rSat ] ; 
     Print [ rSat//N, "     0" ]  } ;

Selection of the equations for solving for the neighbor-effect parameters, q, and conditional probabilities at lattice saturation.

chooseNEandCP::usage = "A routine to determine the values of the neighbor-effect parameters, q, and conditional probabilities when r = rSat." ; 
chooseNEandCP := { 
     If [ n1==n2 && k1 >0 && k2 >0, 
          If [ rSat==1/n1 && w11 >0 && w22 >0, CP:=CP1; CPchoice = "CP1"; 
               If [ w11==w22==w12==w21==1 , (* THEN *) 
                    NE:=NE1a; NEchoice = "NE1a", 
                     (* ELSE *) 
                     If [ !w11==w22, 
                          If [ k1==k2 && w12==w21, NE:=NE1b; NEchoice = "NE1b", 
                               If [ Abs [ deltaA ] > 10^(-12), 
                                    NE:=NE1; NEchoice = "NE1", 
                                    NE:=NE1c; NEchoice = "NE1c" ] ], 
                          If [ Abs [ deltaA ] > 10^(-12), 
                               NE:=NE2a; NEchoice = "NE2a", 
                               NE:=NE2b; NEchoice = "NE2b" ]  ]  ]  ] ; 
          If [ rSat==1/n1 && w11==0 && w22==0, 
               NE:=NE2c; NEchoice = "NE2c"; CP:=CP2; CPchoice = "CP2" ] ; 
          If [ rSat==1/n1 && w11==0 && !w22==0, 
               NE:=NE3; NEchoice = "NE3"; CP:=CP3; CPchoice = "CP3" ] ; 
          If [ rSat==1/n1 && !w11==0 && w22==0, 
               NE:=NE4; NEchoice = "NE4"; CP:=CP4; CPchoice = "CP4" ] ; 
          If [ rSat==1/(n1+1), 
               NE:=NE5; NEchoice = "NE5"; CP:=CP5; CPchoice = "CP5" ] ; 
          If [ rSat==2/(2n1+1) && w22==0, 
               NE:=NE6; NEchoice = "NE6"; CP:=CP6a; CPchoice = "CP6a" ] ; 
          If [ rSat==2/(2n1+1) && w11==0, 
               NE:=NE6; NEchoice = "NE6"; CP:=CP6b; CPchoice = "CP6b" ] ] ; 
     If [ k1==0, NE:=NE7; NEchoice = "NE7"; CP:=CP7; CPchoice = "CP7" ] ; 
     If [ k2==0, NE:=NE8; NEchoice = "NE8"; CP:=CP8; CPchoice = "CP8" ] ; 
     If [ !n1==n2 && k1 >0 && k2 >0, 
          If [ rSat==1/n1 && !( rSat==1/(n2+1) || rSat==2/(n1+n2) || 
               rSat==2/(n1+n2+1) ), 
               NE:=NE9; NEchoice = "NE9"; CP:=CP9; CPchoice = "CP9" ] ; 
          If [ rSat==1/n1 && n1-n2==1 && w11==w22==w12==0, 
               NE:=NE9b; NEchoice = "NE9b"; CP:=CP9b; CPchoice = "CP9b" ] ; 
          If [ rSat==1/n2 && !( rSat==1/(n1+1) || rSat==2/(n1+n2) || 
               rSat==2/(n1+n2+1) ), 
               NE:=NE10; NEchoice = "NE10"; CP:=CP10; CPchoice = "CP10" ] ; 
          If [ rSat==1/n2 && n2-n1==1 && w11==w22==w21==0, 
               NE:=NE10b; NEchoice = "NE10b"; 
               CP:=CP10b; CPchoice = "CP10b" ] ; 
          If [ Abs [ n1-n2 ] == 1 && w11==w12==w21==0 && w22 >0, 
               If [ n1 > n2,  NE:=NE11; NEchoice = "NE11"; 
                    CP:=CP11; CPchoice = "CP11", 
                    NE:=NE12; NEchoice = "NE12"; 
                    CP:=CP12; CPchoice = "CP12" ]  ]  ; 
          If [ Abs [ n1-n2 ] == 1 && w22==w12==w21==0 && w11 >0, 
               If [ n1 > n2,  NE:=NE13; NEchoice = "NE13"; 
                    CP:=CP13; CPchoice = "CP13", 
                    NE:=NE14; NEchoice = "NE14"; 
                    CP:=CP14; CPchoice = "CP14" ]  ]  ; 
          If [ n2-n1==2 && w21==0 && w11 >0 && w22 >0, 
               NE:=NE15; NEchoice = "NE15"; CP:=CP15; CPchoice = "CP15" ] ; 
          If [ n1-n2==2 && w12==0 && w11 >0 && w22 >0, 
               NE:=NE16; NEchoice = "NE16"; CP:=CP16; CPchoice = "CP16" ] ; 
          If [ n1-n2==1 && w12==0 && w11 >0 && w22 >0, 
               NE:=NE17; NEchoice = "NE17"; CP:=CP17; CPchoice = "CP17" ] ; 
          If [ n2-n1==1 && w21==0 && w11 >0 && w22 >0, 
               NE:=NE18; NEchoice = "NE18"; CP:=CP18; CPchoice = "CP18" ] ; 
          If [ n1-n2==2 && w12==0 && (w11==0 || w22==0), 
               NE:=NE19; NEchoice = "NE19" ] ; 
          If [ n2-n1==2 && w21==0 && (w11==0 || w22==0), 
               NE:=NE20; NEchoice = "NE20" ] ; 
          If [ n1-n2 > 2 && w12==0, NE:=NE19; NEchoice = "NE19" ] ; 
          If [ n2-n1 > 2 && w21==0, NE:=NE20; NEchoice = "NE20" ] ; 
          If [ n1-n2 >1 && w11==w12==0 && w22 >0, 
               CP:=CP19; CPchoice = "CP19" ] ; 
          If [ n2-n1 >1 && w11==w21==0 && w22 >0, 
               CP:=CP19; CPchoice = "CP19" ] ; 
          If [ n1-n2 >1 && w22==w12==0 && w11 >0, 
               CP:=CP20; CPchoice = "CP20" ] ; 
          If [ n2-n1 >1 && w22==w21==0 && w11 >0, 
               CP:=CP20; CPchoice = "CP20" ] ; 
          If [ n1-n2 >2 && w12==0 && w11 >0 && w22>0, 
               CP:=CP21; CPchoice = "CP21" ]  ; 
          If [ n2-n1 >2 && w21==0 && w11 >0 && w22>0, 
               CP:=CP21; CPchoice = "CP21" ]  ; 
          If [ n1-n2 >1 && w11==w22==w12==0, 
               CP:=CP22; CPchoice = "CP22" ] ; 
          If [ n2-n1 >1 && w11==w22==w21==0, 
               CP:=CP22; CPchoice = "CP22" ] ; 
          If [ w11==w22==w12==w21==0, CP:=CP22; CPchoice = "CP22"; 
               If [ n1 > n2, NE:=NE19; NEchoice = "NE19", 
                    NE:=NE20; NEchoice = "NE20" ]  ]  ; 
          If [ n1-n2==1 && w11==w12==0 && w22 >0 && w21 >0, 
               NE:=NE21; NEchoice = "NE21"; CP:=CP23; CPchoice = "CP23" ] ; 
          If [ n1-n2==1 && w22==w12==0 && w11 >0 && w21 >0, 
               NE:=NE22; NEchoice = "NE22"; CP:=CP24; CPchoice = "CP24" ] ; 
          If [ n2-n1==1 && w11==w21==0 && w22 >0 && w12 >0, 
               NE:=NE23; NEchoice = "NE23"; CP:=CP25; CPchoice = "CP25" ] ; 
          If [ n2-n1==1 && w22==w21==0 && w11 >0 && w12 >0, 
               NE:=NE24; NEchoice = "NE24"; CP:=CP26; CPchoice = "CP26" ]  ] ; 
     If [ NEchoice=="NE1" || NEchoice=="NE2a" || NEchoice=="NE2c", 
          lastPoint = False, lastPoint = True ] ; 
     r = rSat ; 
     If [ flag, 
          If [ lastPoint, NE ; CP, {E1L,E1R,E2L,E2R,Q} = {1,1,1,1,0.5} ] , 
          NE ; CP ; If [ lastPoint==False, checkLimits ] ]  } ;

Equations for calculating values of neighbor-effect parameters and q when r is at its maximum value; parameter values for special cases

NE1::usage = "Equations for neighbor-effect parameters when r = rSat = 1/n1, both k1 and k2 are > zero, both w11 and w22 are > zero, w11 does not = w22, deltaA does not = 0, and n1 = n2." ;

NE1 := {{E1L, E1R, E2L, E2R, Q} = {e1l, e1r, e2l, e2r, q} /. 
            FindRoot[{
                e1l == ((w11 (1 - q) e2r + 
                      w21 q e1r)/((1 - q) e2r + q e1r)), 
                e1r == ((w11 q e2l + w12 (1 - q) e1l)/(q e2l + 
                      (1 - q) e1l)), 
                e2l == ((w12 (1 - q) e2r + 
                      w22 q e1r)/((1 - q) e2r + q e1r)), 
                e2r == ((w21 q e2l + w22 (1 - q) e1l)/(q e2l + 
                      (1 - q) e1l)), 
                k1 (e1l - e2r) == k2 (e1r - e2l)}, {e1l, E1L}, {
                e1r, E1R}, {e2l, E2L}, {e2r, E2R}, {q, Q}, 
              AccuracyGoal -> 7, MaxIterations -> 200, 
              WorkingPrecision -> 17]};

NE1a::usage = "Equations for neighbor-effect parameters when r = rSat = 1/n1, both k1 and k2 are > zero, w11 = w22 = w12 = w21 = 1, and n1 = n2." ; 
NE1a := {E1L,E1R,E2L,E2R,Q} = {1,1,1,1, k1/(k1 + k2) } ;

NE1b::usage = "Equations for neighbor-effect parameters when r = rSat = 1/n1, k1 = k2, w12 = w21, n1 = n2, both w11 and w22 are > zero, and w11 does not = w22." ;

NE1b := {
        E1L = ((w11 (w12 - w22) + 
              (w11 - w12) (((w11 w22))^0.5))/(w11 - w22)); 
        E2L = ((w22 (w12 - w11) + 
              (w22 - w12) (((w11 w22))^0.5))/(w22 - w11)); {E1R, 
            E2R, Q} = {E1L, E2L, 0.5}};

NE1c::usage = "Equations for neighbor-effect parameters when r = rSat = 1/n1, both k1 and k2 are > zero, both w11 and w22 are > zero, w11 does not = w22, deltaA = 0, and n1 = n2." ;

NE1c := {
        E1L = (((k1 + k2) w11 w)/((k1 + k2) w + 
              k1 (w11 - w21))); 
        E2L = (((k1 + k2) w - k1 E1L)/k2); 
        E1R = (((k1 + k2) w11 w)/((k1 + k2) w + 
              k2 (w11 - w12))); 
        E2R = (((k1 + k2) w - k2 E1R)/k1); 
        Q = ((k1 E1L E2R)/(k1 E1L E2R + k2 E1R E2L))};

NE2a::usage = "Equations for neighbor-effect parameters when r = rSat = 1/n1, both k1 and k2 are > zero, w11 = w22 > 0, w11, w22, w12 and w21 are not all one, deltaA does not = 0, and n1 = n2." ;

NE2a := {{E1L, E1R, E2L, E2R, Q} = {e1l, e1r, e1r, e1l, q} /. 
            FindRoot[{
                e1l == ((w11 (1 - q) e1l + 
                      w21 q e1r)/((1 - q) e1l + q e1r)), 
                e1r == ((w11 q e1r + w12 (1 - q) e1l)/(q e1r + 
                      (1 - q) e1l)), 
                ((k1 (1 - q) e1l)/e1r) == ((k2 q e1r)/e1l)}, {e1l, 
                E1L}, {e1r, E1R}, {q, Q}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 17]};

NE2b::usage = "Equations for neighbor-effect parameters when r = rSat = 1/n1, both k1 and k2 are > zero, w11 = w22 > 0, w11, w22, w12 and w21 are not all one, deltaA = 0, and n1 = n2." ;

NE2b := {Q = (k1/(k1 + k2)); 
        E1L = E1R = 
            E2L = E2R = ((k2 w11 + k1 w21)/(k1 + k2))};

NE2c::usage = "Equations for neighbor-effect parameters when r = rSat = 1/n1, both k1 and k2 are > zero, both w11 and w22 are zero, and n1 = n2." ; 
NE2c := { 
     If [ k1 w21 > k2 w12 , Q = 1 ; E1L = w21 ; E1R = 0 ; 
          E2L = 0 ; E2R = w21 ] ; 
     If [ k1 w21 < k2 w12 , Q = 0 ; E1L = 0 ; E1R = w12 ; 
          E2L = w12 ; E2R = 0 ] ; 
     If [ k1 w21 == k2 w12, Q = k1/(k1 + k2) ; 
          E1L = w21 Q ; E1R = w21 Q ; E2L = w21 Q ; E2R = w21 Q ]  } ;

NE3::usage = "Equations for neighbor-effect parameters when r = rSat = 1/n1, both k1 and k2 are > zero, w11 = zero, w22 is > zero, and n1 = n2." ;

NE3 := {
        If[k1 w21 > k2 w12, Q = 1; E1R = 0; E2R = w21; 
          E1L = w21 - ((w21 w22)/(w22 + ((w21 k1)/k2) - w12)); 
          E2L = ((w21 k1 w22)/(k2 (w22 + ((w21 k1)/k2) - w12)))]; 
        If[k1 w21 < k2 w12, Q = 0; E1L = 0; E2L = w12; 
          E1R = w12 - ((w12 w22)/(w22 + ((w12 k2)/k1) - w21)); 
          E2R = ((w12 k2 w22)/(k1 (w22 + ((w12 k2)/k1) - w21)))]; 
        If[k1 w21 == k2 w12, Q = 0.5; E1L = 0; E1R = 0; E2L = w12; 
          E2R = w21]};

NE4::usage = "Equations for neighbor-effect parameters when r = rSat = 1/n1, both k1 and k2 are > zero, w22 = zero, w11 > zero, and n1 = n2." ;

NE4 := {
        If[k1 w21 > k2 w12, Q = 1; E1L = w21; E2L = 0; 
          E1R = ((w21 k1 w11)/(k2 (w11 + ((w21 k1)/k2) - w12))); 
          E2R = w21 - ((w21 w11)/(w11 + ((w21 k1)/k2) - w12))]; 
        If[k1 w21 < k2 w12, Q = 0; E1R = w12; E2R = 0; 
          E1L = ((w12 k2 w11)/(k1 (w11 + ((w12 k2)/k1) - w21))); 
          E2L = w12 - ((w12 w11)/(w11 + ((w12 k2)/k1) - w21))]; 
        If[k1 w21 == k2 w12, Q = 0.5; E1L = w21; E1R = w12; E2L = 0; 
          E2R = 0]};

NE5::usage = "Neighbor-effect parameters when r = rSat = 1/(n1+1), w11 = w22 = w12 = w21 = zero, both k1 and k2 are > zero, and n1 = n2." ; 
NE5 := {E1L,E1R,E2L,E2R,Q} = {0,0,0,0, k1/(k1+k2) } ;

NE6::usage = "Neighbor-effect parameters when r = rSat = 2/(2n1+1) and n1 = n2." ; 
NE6 := {E1L,E1R,E2L,E2R,Q} = {0,0,0,0,0.5} ;

NE7::usage = "Neighbor-effect parameters when r = rSat and k1 = zero." ; 
NE7 := {E1R,E2L,Q} = {w12,w12,0} ;

NE8::usage = "Neighbor-effect parameters when r = rSat and k2 = zero." ; 
NE8 := {E1L,E2R,Q} = {w21,w21,1} ;

NE9::usage = "Neighbor-effect parameters when r = rSat = 1/n1, k1 and k2 are > zero, n1 < n2, rSat does not = 1/(n2+1), rSat does not = 2/(n1+n2), and rSat does not = 2/(n1+n2+1)." ; 
NE9 := {E1L,E1R,E2L,E2R,Q} = {w21,w11,w22,w21,1} ;

NE9b::usage = "Neighbor-effect parameters when r = rSat = 1/n1 = 1/(n2+1), k1 and k2 are > zero, and w11 = w22 = w12 = 0." ; 
NE9b := { 
     If [ k2/k1==w21, {E1L,E1R,E2L,E2R,Q} = {0,0,0,0,0.5} , 
          If [ k2/k1 > w21, {E1L,E1R,E2L,E2R,Q} = {0,0,0,0,0} , 
               {E1L,E1R,E2L,E2R,Q} := 
                    {w21 (1-FB1R), 0,0, w21 (1-FB1R), 1} ] ] } ;

NE10::usage = "Neighbor-effect parameters when r = rSat = 1/n2, k1 and k2 are > zero, n1 > n2, rSat does not = 1/(n1+1), rSat does not = 2/(n1+n2), and rSat does not = 2/(n1+n2+1)." ; 
NE10 := {E1L,E1R,E2L,E2R,Q} = {w11,w12,w12,w22,0} ;

NE10b::usage = "Neighbor-effect parameters when r = rSat = 1/n2 = 1/(n1+1), k1 and k2 are > zero, and w11 = w22 = w21 = 0." ; 
NE10b := { 
     If [ k1/k2==w12, {E1L,E1R,E2L,E2R,Q} = {0,0,0,0,0.5} , 
          If [ k1/k2 > w12, {E1L,E1R,E2L,E2R,Q} = {0,0,0,0,1} , 
               {E1L,E1R,E2L,E2R,Q} := 
                    {0, w12 (1-FB1L), w12 (1-FB1L), 0, 0} ] ] } ;

NE11::usage = "Neighbor-effect parameters when r = rSat = 1/(n2 + 1) = 2/(n1 + n2 + 1), n1 - n2 = 1, k1 and k2 are > zero, w11 = w12 = w21 = zero, and w22 > zero." ;

NE11 := {
        Q = 1 - (1/2) (1 + (1/(((1 + ((4 w22 k1)/k2)))^0.5))); 
        E2R = ((w22 (1 - 2 Q))/(1 - Q)); {E1L, E1R, E2L} = {0, 0, 
            0}};

NE12::usage = "Neighbor-effect parameters when r = rSat = 1/(n1 + 1) = 2/(n1 + n2 + 1), n2 - n1 = 1, k1 and k2 are > zero, w11 = w12 = w21 = zero, and w22 > zero." ;

NE12 := {Q = (1/2) (1 + (1/(((1 + ((4 w22 k2)/k1)))^0.5))); 
        E2L = ((w22 (2 Q - 1))/Q); {E1L, E1R, E2R} = {0, 0, 0}};

NE13::usage = "Neighbor-effect parameters when r = rSat = 1/(n2 + 1) = 2/(n1 + n2 + 1), n1 - n2 = 1, k1 and k2 are > zero, w22 = w12 = w21 = zero, and w11 > zero." ;

NE13 := {
        Q = 1 - (1/2) (1 + (1/(((1 + ((4 w11 k1)/k2)))^0.5))); 
        E1L = ((w11 (1 - 2 Q))/(1 - Q)); {E1R, E2L, E2R} = {0, 0, 
            0}};

NE14::usage = "Neighbor-effect parameters when r = rSat = 1/(n1 + 1) = 2/(n1 + n2 + 1), n2 - n1 = 1, k1 and k2 are > zero, w22 = w12 = w21 = zero, and w11 > zero." ;

NE14 := {Q = (1/2) (1 + (1/(((1 + ((4 w11 k2)/k1)))^0.5))); 
        E1R = ((w11 (2 Q - 1))/Q); {E1L, E2L, E2R} = {0, 0, 0}};

NE15::usage = "Neighbor-effect parameters when r = rSat = 1/(n1 + 1) = 2/(n1 + n2), n2 - n1 = 2, k1 and k2 are > zero, w21 = zero, and both w11 and w22 are > zero." ;

NE15 := {c = (Sqrt[((4 w11 w22 k2)/k1) + 1]); 
        Q = ((c + 1)/(2 c)); {E1L, E1R, E2L, E2R} = {0, w11, w22, 0}};

NE16::usage = "Neighbor-effect parameters when r = rSat = 1/(n2 + 1) = 2/(n1 + n2), n1 - n2 = 2, k1 and k2 are > zero, w12 = zero, and both w11 and w22 are > zero." ;

NE16 := {c = (Sqrt[((4 w11 w22 k1)/k2) + 1]); 
        Q = ((c - 1)/(2 c)); {E1L, E1R, E2L, E2R} = {w11, 0, 0, w22}};

NE17::usage = "Neighbor-effect parameters when r = rSat = 2/(n1 + n2), n1 - n2 = 1, k1 and k2 are > zero, w12 = zero, and both w11 and w22 are > zero." ; 
NE17 := {E1L,E1R,E2L,E2R,Q} = {w11,0,0,w22,0.5} ;

NE18::usage = "Neighbor-effect parameters when r = rSat = 2/(n1 + n2), n2 - n1 = 1, k1 and k2 are > zero, w21 = zero, and both w11 and w22 are > zero." ; 
NE18 := {E1L,E1R,E2L,E2R,Q} = {0,w11,w22,0,0.5} ;

NE19::usage = "Neighbor-effect parameters when r = rSat = 1/(n2 + 1), k1 and k2 are > zero, and one of the following is true: 
a)  n1- n2 = 1 and w11, w22, w12, and w21 are all zero;
b)  n1- n2 = 2, w12 = zero and at least one of w11 or w22 is zero;
c)  n1- n2 >2 and w12 is zero." ; 
NE19 := {E1L,E1R,E2L,E2R,Q} = {w11,0,0,w22,0} ;

NE20::usage = "Neighbor-effect parameters when r = rSat = 1/(n1 + 1), k1 and k2 are > zero, and one of the following is true: 
a)  n2- n1 = 1 and w11, w22, w12, and w21 are all zero;
b)  n2- n1 = 2, w21 = zero and at least one of w11 or w22 is zero;
c)  n2- n1 >2 and w21 is zero." ; 
NE20 := {E1L,E1R,E2L,E2R,Q} = {0,w11,w22,0,1} ;

NE21::usage = "Neighbor-effect parameters when r = rSat = 1/n1 = 1/(n2 + 1) = 2/(n1 + n2 + 1), n1- n2 = 1, k1 and k2 are > zero, w11 = w12 = zero, and w22 and w21 are > zero." ;

NE21 := 
      Block[{b, k, y, z}, k = (k2/k1); 
        z = (1/2) (k - w21 + (Sqrt[(((k - w21))^2) + 4 w22 k])); 
        b = w22 w21 - (z^2); 
        E2R = ((b + 
              (Sqrt[(b^2) + 4 w22 z (z + w22) (z + w21)]))/(2 
              (z + w22))); y = (z/E2R); 
        Q = 1 - (y/(1 + y)); {E1L, E1R, E2L} = {0, 0, 0}];

NE22::usage = "Neighbor-effect parameters when r = rSat = 1/n1 = 1/(n2 + 1) = 2/(n1 + n2 + 1), n1- n2 = 1, k1 and k2 are > zero, w22 = w12 = zero, and w11 and w21 are > zero." ;

NE22 := 
      Block[{b, k, y, z}, k = (k2/k1); 
        z = (1/2) (k - w21 + (Sqrt[(((k - w21))^2) + 4 w11 k])); 
        b = w11 w21 - (z^2); 
        E1L = ((b + 
              (Sqrt[(b^2) + 4 w11 z (z + w11) (z + w21)]))/(2 
              (z + w11))); y = (z/E1L); 
        Q = 1 - (y/(1 + y)); {E1R, E2R, E2L} = {0, 0, 0}];

NE23::usage = "Neighbor-effect parameters when r = rSat = 1/n2 = 1/(n1 + 1) = 2/(n1 + n2 + 1), n2 - n1 = 1, k1 and k2 are > zero, w11 = w21 = zero, and w22 and w12 are > zero." ;

NE23 := 
      Block[{b, k, y, z}, k = (k1/k2); 
        z = (1/2) (k - w12 + (Sqrt[(((k - w12))^2) + 4 w22 k])); 
        b = w22 w12 - (z^2); 
        E2L = ((b + 
              (Sqrt[(b^2) + 4 w22 z (z + w22) (z + w12)]))/(2 
              (z + w22))); y = (z/E2L); 
        Q = (y/(1 + y)); {E1L, E1R, E2R} = {0, 0, 0}];

NE24::usage = "Neighbor-effect parameters when r = rSat = 1/n2 = 1/(n1 + 1) = 2/(n1 + n2 + 1), n2 - n1 = 1, k1 and k2 are > zero, w22 = w21 = zero, and w11 and w12 are > zero." ;

NE24 := 
      Block[{b, k, y, z}, k = (k1/k2); 
        z = (1/2) (k - w12 + (Sqrt[(((k - w12))^2) + 4 w11 k])); 
        b = w11 w12 - (z^2); 
        E1R = ((b + 
              (Sqrt[(b^2) + 4 w11 z (z + w11) (z + w12)]))/(2 
              (z + w11))); y = (z/E1R); 
        Q = (y/(1 + y)); {E1L, E2L, E2R} = {0, 0, 0}];

Expressions for conditional probabilities when r 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." ; 
b1Lf::usage = "b1Lf is the probability that a randomly chosen b1L end of a bound ligand is bordered by a free lattice residue." ; 
b1Rf::usage = "b1Rf is the probability that a randomly chosen b1R end of a bound ligand is bordered by a free lattice residue." ; 
b2Lf::usage = "b2Lf is the probability that a randomly chosen b2L end of a bound ligand is bordered by a free lattice residue." ; 
b2Rf::usage = "b2Rf is the probability that a randomly chosen b2R end of a bound ligand is bordered by a free lattice residue." ; 
fb1L::usage = "fb1L is the probability that a randomly chosen free lattice residue is bordered on the right by a b1L end of a bound ligand." ; 
fb1R::usage = "fb1R is the probability that a randomly chosen free lattice residue is bordered on the left by a b1R end of a bound ligand." ; 
fb2L::usage = "fb2L is the probability that a randomly chosen free lattice residue is bordered on the right by a b2L end of a bound ligand." ; 
fb2R::usage = "fb2R is the probability that a randomly chosen free lattice residue is bordered on the left by a b2R end of a bound ligand." ; 
b1Lb1R::usage = "b1Lb1R is the probability that a randomly chosen b1L end of a bound ligand is bordered by a b1R end of a bound ligand." ; 
b1Rb1L::usage = "b1Rb1L is the probability that a randomly chosen b1R end of a bound ligand is bordered by a b1L end of a bound ligand." ; 
b2Lb2R::usage = "b2Lb2R is the probability that a randomly chosen b2L end of a bound ligand is bordered by a b2R end of a bound ligand." ; 
b2Rb2L::usage = "b2Rb2L is the probability that a randomly chosen b2R end of a bound ligand is bordered by a b2L end of a bound ligand." ; 
b1Lb2R::usage = "b1Lb2R (= b2Rb1L) is the probability that a randomly chosen b1L end of a bound ligand is bordered by a b2R end of a bound ligand (or vice versa)." ; 
b1Rb2L::usage = "b1Rb2L (= b2Lb1R) is the probability that a randomly chosen b1R end of a bound ligand is bordered by a b2L end of a bound ligand (or vice versa)." ;

ff := 1/(1 + Q X/E1L + (1-Q) X/E2L ) ; 
b1Lf [h_] := h[[8]]/h[[2]] /;  h[[1]] > 0 ; (* = ff/E1L *) 
b1Lf [h_] := B1LF /;  h[[1]]==0 ; 
b1Rf [h_] := h[[8]]/h[[3]] /;  h[[1]] > 0 ; (* = ff/E1R *) 
b1Rf [h_] := B1RF /;  h[[1]]==0 ; 
b2Lf [h_] := h[[8]]/h[[4]] /;  h[[1]] > 0 ; (* = ff/E2L *) 
b2Lf [h_] := B2LF /;  h[[1]]==0 ; 
b2Rf [h_] := h[[8]]/h[[5]] /;  h[[1]] > 0 ; (* = ff/E2R *) 
b2Rf [h_] := B2RF /;  h[[1]]==0 ; 
fb1L [h_] := h[[6]] h[[7]] h[[8]]/h[[2]] /;  h[[1]] > 0 ; (* = Q X ff/E1L *) 
fb1L [h_] := FB1L /;  h[[1]]==0 ; 
fb1R [h_] := (1-h[[6]]) h[[7]] h[[8]]/h[[3]] /;  h[[1]] > 0 ; (* = (1-Q) X ff/E1R *) 
fb1R [h_] := FB1R /;  h[[1]]==0 ; 
fb2L [h_] := (1-h[[6]]) h[[7]] h[[8]]/h[[4]] /;  h[[1]] > 0 ; (* = (1-Q) X ff/E2L *) 
fb2L [h_] := FB2L /;  h[[1]]==0 ; 
fb2R [h_] := h[[6]] h[[7]] h[[8]]/h[[5]] /;  h[[1]] > 0 ; (* = Q X ff/E2R *) 
fb2R [h_] := FB2R /;  h[[1]]==0 ; 
b1Lb1R [h_] := (1-h[[6]]) h[[7]] w11 h[[8]]/(h[[2]] h[[3]]) /;  h[[1]] > 0 ; 
     (* = (1-Q) X w11 ff/(E1L E1R) *) 
b1Lb1R [h_] := B1LB1R /;  h[[1]]==0 ; 
b1Rb1L [h_] := h[[6]] h[[7]] w11 h[[8]]/(h[[3]] h[[2]]) /;  h[[1]] > 0 ; 
     (* = Q X w11 ff/(E1R E1L) *) 
b1Rb1L [h_] := B1RB1L /;  h[[1]]==0 ; 
b2Lb2R [h_] := h[[6]] h[[7]] w22 h[[8]]/(h[[4]] h[[5]]) /;  h[[1]] > 0 ; 
     (* = Q X w22 ff/(E2L E2R) *) 
b2Lb2R [h_] := B2LB2R /;  h[[1]]==0 ; 
b2Rb2L [h_] := (1-h[[6]]) h[[7]] w22 h[[8]]/(h[[5]] h[[4]]) /;  h[[1]] > 0 ; 
     (* = (1-Q) X w22 ff/(E2R E2L) *) 
b2Rb2L [h_] := B2RB2L /;  h[[1]]==0 ; 
b1Lb2R [h_] := h[[6]] h[[7]] w21 h[[8]]/(h[[2]] h[[5]]) /;  h[[1]] > 0 ; 
     (* = Q X w21 ff/(E1L E2R) *) 
b1Lb2R [h_] := B1LB2R /;  h[[1]]==0 ; 
b1Rb2L [h_] := (1-h[[6]]) h[[7]] w12 h[[8]]/(h[[3]] h[[4]]) /;  h[[1]] > 0 ; 
     (* = (1-Q) X w12 ff/(E1R E2L) *) 
b1Rb2L [h_] := B1RB2L /;  h[[1]]==0 ;

Expressions for conditional probabilities when r is at its maximum value;  probabilities for special cases

B1LF::usage = "B1LF is the value of b1Lf when the lattice is saturated." ; 
B1RF::usage = "B1RF is the value of b1Rf when the lattice is saturated." ; 
B2LF::usage = "B2LF is the value of b2Lf when the lattice is saturated." ; 
B2RF::usage = "B2RF is the value of b2Rf when the lattice is saturated." ; 
FB1L::usage = "FB1L is the value of fb1L when the lattice is saturated." ; 
FB1R::usage = "FB1R is the value of fb1R when the lattice is saturated." ; 
FB2L::usage = "FB2L is the value of fb2L when the lattice is saturated." ; 
FB2R::usage = "FB2R is the value of fb2R when the lattice is saturated." ; 
B1LB1R::usage = "B1LB1R is the value of b1Lb1R when the lattice is saturated." ; 
B1RB1L::usage = "B1RB1L is the value of b1Rb1L when the lattice is saturated." ; 
B2LB2R::usage = "B2LB2R is the value of b2Lb2R when the lattice is saturated." ; 
B2RB2L::usage = "B2RB2L is the value of b2Rb2L when the lattice is saturated." ; 
B1LB2R::usage = "B1LB2R is the value of b1Lb2R when the lattice is saturated." ; 
B1RB2L::usage = "B1RB2L is the value of b1Rb2L when the lattice is saturated." ;

CP1::usage = "Expressions for conditional probabilities when r =  rSat = 1/n1, both k1 and k2 are > zero, and n1 = n2." ;

CP1 := {B1LF = 0; B1RF = 0; B2LF = 0; B2RF = 0; 
        FB1L = (Q/(((Q/E1L) + ((1 - Q)/E2L)) E1L)); 
        FB1R = ((1 - Q)/(((Q/E1L) + ((1 - Q)/E2L)) E1R)); 
        FB2L = ((1 - Q)/(((Q/E1L) + ((1 - Q)/E2L)) E2L)); 
        FB2R = (Q/(((Q/E1L) + ((1 - Q)/E2L)) E2R)); 
        B1LB1R = 
          (((1 - Q) w11)/(((Q/E1L) + ((1 - Q)/E2L)) E1L E1R)); 
        B1RB1L = ((Q w11)/(((Q/E1L) + ((1 - Q)/E2L)) E1R E1L)); 
        B2LB2R = ((Q w22)/(((Q/E1L) + ((1 - Q)/E2L)) E2L E2R)); 
        B2RB2L = 
          (((1 - Q) w22)/(((Q/E1L) + ((1 - Q)/E2L)) E2R E2L)); 
        B1LB2R = ((Q w21)/(((Q/E1L) + ((1 - Q)/E2L)) E1L E2R)); 
        B1RB2L = 
          (((1 - Q) w12)/(((Q/E1L) + ((1 - Q)/E2L)) E1R E2L))};

CP2::usage = "Values of conditional probabilities when r = rSat = 1/n1, both k1 and k2 are > zero, w11 and w22 are zero, and n1 = n2." ; 
CP2 := { 
     If [ k1 w21 > k2 w12 , 
          B1RB2L = k2 w12/(k1 w21) ; B1RF = 1 - B1RB2L ; 
          B2LF = 1 - B1RB2L ; B1LF = 0 ; B2RF = 0 ; FB1L = 1 ; 
          FB1R = 0 ; FB2L = 0 ; FB2R = 1 ; B1LB1R = 0 ; 
          B1RB1L = 0 ; B2LB2R = 0 ; B1LB2R = 1 ; B2RB2L = 0 ] ; 
     If [ k1 w21 < k2 w12 , 
          B1LB2R = k1 w21/(k2 w12) ; B1LF = 1 - B1LB2R ; 
          B2RF = 1 - B1LB2R ; B1RF = 0 ; B2LF = 0 ; FB1L = 0 ; 
          FB1R = 1 ; FB2L = 1 ; FB2R = 0 ; B1LB1R = 0 ; 
          B1RB1L = 0 ; B2LB2R = 0 ; B1RB2L = 1 ; B2RB2L = 0 ] ; 
     If [ k1 w21 == k2 w12, CP1 ]  } ;

CP3::usage = "Values of conditional probabilities when r = rSat = 1/n1, both k1 and k2 are > zero, w11 = zero, w22 > zero, and n1 = n2." ; 
CP3 := { 
     B2LF = 0 ; B2RF = 0 ; 
     If [ k1 w21 > k2 w12 , 
          B2LB2R = (w22 E1L)/(E2L E2R) ; 
          B1RF = B2LB2R ; FB2R = E1L/E2R ; FB1R = 1- FB2R ; 
          B1LF = 0 ; FB1L = 1 ; FB2L = 0 ; B1LB1R = 0 ; 
          B1RB1L = 0 ; B2RB2L = 0 ; B1LB2R = 1 ; 
          If [ w12==0, B1RB2L = 0, B1RB2L = 1 - B2LB2R ]  ] ; 
     If [ k1 w21 < k2 w12 , 
          B2RB2L = (w22 E1R)/(E2R E2L) ; 
          B1LF = B2RB2L ; FB2L = E1R/E2L ; FB1L = 1- FB2L ; 
          B1RF = 0 ; FB1R = 1 ; FB2R = 0 ; B1RB1L = 0 ; 
          B1LB1R = 0 ; B2LB2R = 0 ; B1RB2L = 1 ; 
          If [ w21==0, B1LB2R = 0, B1LB2R = 1 - B2RB2L ]  ] ; 
     If [ k1 w21 == k2 w12, 
          B1LF = 0 ; B1RF = 0 ; FB1L = 1 ; FB1R = 1 ; FB2L = 0 ; 
          FB2R = 0 ; B1LB1R = 0 ; B1RB1L = 0 ; B2LB2R = 0 ; 
          B2RB2L = 0 ; B1RB2L = 1 ; B1LB2R = 1 ]  }  ;

CP4::usage = "Values of conditional probabilities when r = rSat = 1/n1, both k1 and k2 are > zero, w22 = zero, w11 > zero, and n1 = n2." ; 
CP4 := {  
     B1LF = 0 ; B1RF = 0 ; 
     If [ k1 w21 > k2 w12 , 
          B1RB1L = (w11 E2R)/(E1R E1L) ; 
          B2LF = B1RB1L ; FB1L = E2R/E1L ; FB2L = 1- FB1L ; 
          B2RF = 0 ; FB1R = 0 ; FB2R = 1 ; B1LB1R = 0 ; 
          B2RB2L = 0 ; B2LB2R = 0 ; B1LB2R = 1 ; 
          If [ w12==0, B1RB2L = 0, B1RB2L = 1 - B1RB1L ]  ] ; 
     If [ k1 w21 < k2 w12 , 
          B1LB1R = (w11 E2L)/(E1R E1L) ; 
          B2RF = B1LB1R ; FB1R = E2L/E1R ; FB2R = 1- FB1R ; 
          B2LF = 0 ; FB1L = 0 ; FB2L = 1 ; B1RB1L = 0 ; 
          B2LB2R = 0 ; B2RB2L = 0 ; B1RB2L = 1 ; 
          If [ w21==0, B1LB2R = 0, B1LB2R = 1 - B1LB1R ]  ] ; 
     If [ k1 w21 == k2 w12, 
          B2LF = 0 ; B2RF = 0 ; FB1L = 0 ; FB1R = 0 ; 
          FB2L = 1 ; FB2R = 1 ; B1LB1R = 0 ; B1RB1L = 0 ; 
          B2LB2R = 0 ; B2RB2L = 0 ; B1RB2L = 1 ; B1LB2R = 1 ]  } ;

CP5::usage = "Values of conditional probabilities when r =  rSat = 1/(n1+1), w11 = w22 = w12 = w21 = zero, both k1 and k2 are > zero, and n1 = n2." ; 
CP5 := {
     B1LF = 1 ; B1RF = 1 ; B2LF = 1 ; B2RF = 1 ; 
     FB1L = Q ; FB1R = 1-Q ; FB2L = 1-Q ; FB2R = Q ; 
     B1LB1R = 0 ; B1RB1L = 0 ; B2LB2R = 0 ; B2RB2L = 0 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP6a::usage = "Values of conditional probabilities when r = rSat = 2/(2n1+1), w22 = w12 = w21 = zero, and n1 = n2." ; 
CP6a := { 
     B1LF = 0 ; B1RF = 0 ; B2LF = 1 ; B2RF = 1 ; 
     FB1L = 0 ; FB1R = 0 ; FB2L = 1 ; FB2R = 1 ; 
     B1LB1R = 1 ; B1RB1L = 1 ; B2LB2R = 0 ; B2RB2L = 0 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP6b::usage = "Values of conditional probabilities when r = rSat = 2/(2n1+1), w11 = w12 = w21 = zero, and n1 = n2." ; 
CP6b := { 
     B1LF = 1 ; B1RF = 1 ; B2LF = 0 ; B2RF = 0 ; 
     FB1L = 1 ; FB1R = 1 ; FB2L = 0 ; FB2R = 0 ; 
     B1LB1R = 0 ; B1RB1L = 0 ; B2LB2R = 1 ; B2RB2L = 1 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP7::usage = "Values of conditional probabilities when r = rSat and k1 = zero." ; 
CP7 := {
     FB1L = 0 ; FB1R = 1 ; FB2L = 1 ; FB2R = 0 ; 
     B1RB1L = 0 ; B2LB2R = 0 ; 
     If [ w12 == 0, B1RF = 1 ; B2LF = 1 ; B1RB2L = 0,
          B1RF = 0 ; B2LF = 0 ; B1RB2L = 1 ]  } ;

CP8::usage = "Values of conditional probabilities when r = rSat and k2 = zero." ; 
CP8 := {
     FB1L = 1 ; FB1R = 0 ; FB2L = 0 ; FB2R = 1 ; 
     B1LB1R = 0 ; B2RB2L = 0 ; 
     If [ w21 == 0, B1LF = 1 ; B2RF = 1 ; B1LB2R = 0,
          B1LF = 0 ; B2RF = 0 ; B1LB2R = 1 ]  } ;

CP9::usage = "Values of conditional probabilities when r = rSat = 1/n1, k1 and k2 are > zero, n1 < n2, rSat does not = 1/(n2+1), rSat does not = 2/(n1+n2), and rSat does not = 2/(n1+n2+1)." ; 
CP9 := { 
     B1LF = 0 ; B2RF = 0 ; 
     FB1L = 1 ; FB1R = 0 ; FB2L = 0 ; FB2R = 1 ; 
     B1LB1R = 0 ; B2RB2L = 0 ; 
     B1LB2R = 1 ; B1RB2L = 0 ; 
     If [ w11==0, B1RF = 1 ; B1RB1L = 0, B1RF = 0 ; B1RB1L = 1 ] ; 
     If [ w22==0, B2LF = 1 ; B2LB2R = 0, B2LF = 0 ; B2LB2R = 1 ] } ;

CP9b::usage = "Values of conditional probabilities when r = rSat = 1/n1 = 1/(n2+1), k1 and k2 are > zero, and w11 = w22 = w12 = 0." ; 
CP9b := { 
     B1RF = 1 ; B2LF = 1 ; B1LB1R = 0 ; B1RB1L = 0 ; 
     B2LB2R = 0 ; B2RB2L = 0 ; B1RB2L = 0 ; 
     If [ k2/k1==w21, 
          B1LF = 0 ; B2RF = 0 ; FB1L = 0 ; FB1R = 1 ; 
          FB2L = 1 ; FB2R = 0 ; B1LB2R = 1 , 
          If [ k2/k1 > w21, B1LF := 1-B1LB2R ; B2RF := 1-B1LB2R ; 
               FB1L = 0 ; FB1R = 1 ; FB2L = 1 ; FB2R = 0 ; 
               B1LB2R := w21 k1/k2 , 
               B1LF = 0 ; B2RF = 0 ; FB1L := 1-FB1R ; FB1R = k2/(k1 w21) ; 
               FB2L := FB1R ; FB2R := 1-FB1R ; B1LB2R = 1 ] ] } ;

CP10::usage = "Values of conditional probabilities when r = rSat = 1/n2, k1 and k2 are > zero, n1 < n2, rSat does not = 1/(n1+1), rSat does not = 2/(n1+n2), and rSat does not = 2/(n1+n2+1)." ; 
CP10 := { 
     B1RF = 0 ; B2LF = 0 ; 
     FB1L = 0 ; FB1R = 1 ; FB2L = 1 ; FB2R = 0 ; 
     B1RB1L = 0 ; B2LB2R = 0 ; 
     B1LB2R = 0 ; B1RB2L = 1 ; 
     If [ w11==0, B1LF = 1 ; B1LB1R = 0, B1LF = 0 ; B1LB1R = 1 ] ; 
     If [ w22==0, B2RF = 1 ; B2RB2L = 0, B2RF = 0 ; B2RB2L = 1 ]  } ;

CP10b::usage = "Values of conditional probabilities when r = rSat = 1/n2 = 1/(n1+1), k1 and k2 are > zero, and w11 = w22 = w21 = 0." ; 
CP10b := { 
     B2RF = 1 ; B1LF = 1 ; B2LB2R = 0 ; B2RB2L = 0 ; 
     B1LB1R = 0 ; B1RB1L = 0 ; B1LB2R = 0 ; 
     If [ k1/k2==w12, 
          B2LF = 0 ; B1RF = 0 ; FB2L = 0 ; FB2R = 1 ; 
          FB1L = 1 ; FB1R = 0 ; B1RB2L = 1 , 
          If [ k1/k2 > w12, B2LF := 1-B1RB2L ; B1RF := 1-B1RB2L ; 
               FB2L = 0 ; FB2R = 1 ; FB1L = 1 ; FB1R = 0 ; 
               B1RB2L = w12 k2/k1 , 
               B2LF = 0 ; B1RF = 0 ; FB2L := 1-FB1L ; FB2R := FB1L ; 
               FB1L = k1/(k2 w12) ; FB1R := 1-FB1L ; B1RB2L = 1 ] ] } ;

CP11::usage = "Conditional probabilities when r = rSat = 1/(n2 + 1) = 2/(n1 + n2 + 1), n1 - n2 = 1, k1 and k2 are > zero, w11 = w12 = w21 = zero, and w22 > zero." ; 
CP11 := { 
     B1LF = 1 ; B1RF = 1 ; B2LF = 1-Q/(1-Q) ; B2RF = 0 ; 
     FB1L = Q/(1-Q) ; FB1R = 1 ; FB2L = 1-Q/(1-Q) ; FB2R = 0 ; 
     B1LB1R = 0 ; B1RB1L = 0 ; B2LB2R = Q/(1-Q) ; B2RB2L = 1 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP12::usage = "Conditional probabilities when r = rSat = 1/(n1 + 1) = 2/(n1 + n2 + 1), n2 - n1 = 1, k1 and k2 are > zero, w11 = w12 = w21 = zero, and w22 > zero." ; 
CP12 := { 
     B1LF = 1 ; B1RF = 1 ; B2LF = 0 ; B2RF = 1-(1-Q)/Q ; 
     FB1L = 1 ; FB1R = (1-Q)/Q ; FB2L = 0 ; FB2R = 1-(1-Q)/Q ; 
     B1LB1R = 0 ; B1RB1L = 0 ; B2LB2R = 1 ; B2RB2L = (1-Q)/Q ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP13::usage = "Conditional probabilities when r = rSat = 1/(n2 + 1) = 2/(n1 + n2 + 1), n1 - n2 = 1, k1 and k2 are > zero, w22 = w12 = w21 = zero, and w11 > zero." ; 
CP13 := { 
     B1LF = 0 ; B1RF = 1-Q/(1-Q) ; B2LF = 1 ; B2RF = 1 ; 
     FB1L = 0 ; FB1R = 1-Q/(1-Q) ; FB2L = 1 ; FB2R = Q/(1-Q) ; 
     B1LB1R = 1 ; B1RB1L = Q/(1-Q) ; B2LB2R = 0 ; B2RB2L = 0 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP14::usage = "Conditional probabilities when r = rSat = 1/(n1 + 1) = 2/(n1 + n2 + 1), n2 - n1 = 1, k1 and k2 are > zero, w22 = w12 = w21 = zero, and w11 > zero." ; 
CP14 := { 
     B1LF = 1-(1-Q)/Q ; B1RF = 0 ; B2LF = 1 ; B2RF = 1 ; 
     FB1L = 1-(1-Q)/Q ; FB1R = 0 ; FB2L = (1-Q)/Q ; FB2R = 1 ; 
     B1LB1R = (1-Q)/Q ; B1RB1L = 1 ; B2LB2R = 0 ; B2RB2L = 0 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP15::usage = "Conditional probabilities when r = rSat = 1/(n1 + 1) = 2/(n1 + n2), n2 - n1 = 2, k1 and k2 are > zero, w21 = zero, and both w11 and w22 are > zero." ; 
CP15 := { 
     B1LF = 1-(1-Q)/Q ; B1RF = 0 ; B2LF = 0 ; B2RF = 1-(1-Q)/Q ; 
     FB1L = 1 ; FB1R = 0 ; FB2L = 0 ; FB2R = 1 ; 
     B1LB1R = (1-Q)/Q ; B1RB1L = 1 ; B2LB2R = 1 ; B2RB2L = (1-Q)/Q ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP16::usage = "Conditional probabilities when r = rSat = 1/(n2 + 1) = 2/(n1 + n2), n1 - n2 = 2, k1 and k2 are > zero, w12 = zero, and both w11 and w22 are > zero." ; 
CP16 := { 
     B1LF = 0 ; B1RF = 1-Q/(1-Q) ; B2LF = 1-Q/(1-Q) ; B2RF = 0 ; 
     FB1L = 0 ; FB1R = 1 ; FB2L = 1 ; FB2R = 0 ; 
     B1LB1R = 1 ; B1RB1L = Q/(1-Q) ; B2LB2R = Q/(1-Q) ; B2RB2L = 1 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP17::usage = "Conditional probabilities when r = rSat = 2/(n1 + n2), n1 - n2 = 1, k1 and k2 are > zero, w12 = zero, and both w11 and w22 are > zero." ; 
CP17 := { 
     B1LF = 0 ; B1RF = 0 ; B2LF = 0 ; B2RF = 0 ; 
     FB1L = 0 ; FB1R = 1 ; FB2L = 1 ; FB2R = 0 ; 
     B1LB1R = 1 ; B1RB1L = 1 ; B2LB2R = 1 ; B2RB2L = 1 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP18::usage = "Conditional probabilities when r = rSat = 2/(n1 + n2), n2 - n1 = 1, k1 and k2 are > zero, w21 = zero, and both w11 and w22 are > zero." ; 
CP18 := { 
     B1LF = 0 ; B1RF = 0 ; B2LF = 0 ; B2RF = 0 ; 
     FB1L = 1 ; FB1R = 0 ; FB2L = 0 ; FB2R = 1 ; 
     B1LB1R = 1 ; B1RB1L = 1 ; B2LB2R = 1 ; B2RB2L = 1 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP19::usage = "Conditional probabilities when k1 and k2 are > zero, w11 = zero, w22 > zero, and one of the following is true: 
a)  r = rSat = 1/(n2 + 1), n1- n2 >1, and w12 = zero;
b)  r = rSat = 1/(n1 + 1), n2- n1 >1, and w21 = zero." ; 
CP19 := { 
     B1LF = 1 ; B1RF = 1 ; B2LF = 1-Q ; B2RF = Q ; 
     FB1L = Q ; FB1R = 1-Q ; FB2L = 1-Q ; FB2R = Q ; 
     B1LB1R = 0 ; B1RB1L = 0 ; B2LB2R = Q ; B2RB2L = 1-Q ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP20::usage = "Conditional probabilities when k1 and k2 are > zero, w22 = zero, w11 > zero, and one of the following is true: 
a)  r = rSat = 1/(n2 + 1), n1- n2 >1, and w12 = zero;
b)  r = rSat = 1/(n1 + 1), n2- n1 >1, and w21 = zero." ; 
CP20 := { 
     B1LF = Q ; B1RF = 1-Q ; B2LF = 1 ; B2RF = 1 ; 
     FB1L = Q ; FB1R = 1-Q ; FB2L = 1-Q ; FB2R = Q ; 
     B1LB1R = 1-Q ; B1RB1L = Q ; B2LB2R = 0 ; B2RB2L = 0 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP21::usage = "Conditional probabilities when k1 and k2 are > zero, w11 and w22 are both > zero, and one of the following is true: 
a)  r = rSat = 1/(n2 + 1), n1- n2 >2, and w12 = zero;
b)  r = rSat = 1/(n1 + 1), n2- n1 >2, and w21 = zero." ; 
CP21 := { 
     B1LF = Q ; B1RF = 1-Q ; B2LF = 1-Q ; B2RF = Q ; 
     FB1L = Q ; FB1R = 1-Q ; FB2L = 1-Q ; FB2R = Q ; 
     B1LB1R = 1-Q ; B1RB1L = Q ; B2LB2R = Q ; B2RB2L = 1-Q ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP22::usage = "Conditional probabilities when k1 and k2 are > zero, w11 = w22 = zero, and one of the following is true: 
a)  r = rSat = 1/(n2 + 1), n1- n2 >1, and w12 = zero;
b)  r = rSat = 1/(n1 + 1), n2- n1 >1, and w21 = zero;
c)  r = rSat = 1/(n2 + 1), n1- n2 >0, and w12 = w21 = zero;
d)  r = rSat = 1/(n1 + 1), n2- n1 >0, and w12 = w21 = zero." ; 
CP22 := { 
     B1LF = 1 ; B1RF = 1 ; B2LF = 1 ; B2RF = 1 ; 
     FB1L = Q ; FB1R = 1-Q ; FB2L = 1-Q ; FB2R = Q ; 
     B1LB1R = 0 ; B1RB1L = 0 ; B2LB2R = 0 ; B2RB2L = 0 ; 
     B1LB2R = 0 ; B1RB2L = 0  } ;

CP23::usage = "Conditional probabilities when r = rSat = 1/n1 = 1/(n2 + 1) = 2/(n1 + n2 + 1), n1- n2 = 1, k1 and k2 are > zero, w11 = w12 = zero, and w22 and w21 are > zero." ; 
CP23 := { 
     FB1L = (k1 w21 - k2)/(k1 w21 - k2 (1-Q)/Q) ;  
     B1LF = FB1L (1-Q)/Q ; B1RF = 1 ; B2LF = 1- FB1L ; B2RF = 0 ; 
     FB1R = 1 ; FB2L = 1- FB1L ; FB2R = 0 ; 
     B1LB1R = 0 ; B1RB1L = 0 ; B2LB2R = FB1L ; B2RB2L = B1LF ; 
     B1LB2R = 1- B1LF ; B1RB2L = 0  } ;

CP24::usage = "Conditional probabilities when r = rSat = 1/n1 = 1/(n2 + 1) = 2/(n1 + n2 + 1), n1- n2 = 1, k1 and k2 are > zero, w22 = w12 = zero, and w11 and w21 are > zero." ; 
CP24 := { 
     FB2R = (k1 w21 - k2)/(k1 w21 - k2 (1-Q)/Q) ;  
     B1LF = 0 ; B1RF = 1- FB2R ; B2LF = 1 ; B2RF = FB2R (1-Q)/Q ; 
     FB1L = 0 ; FB1R = 1- FB2R ; FB2L = 1 ; 
     B1LB1R = B2RF ; B1RB1L = FB2R ; B2LB2R = 0 ; B2RB2L = 0 ; 
     B1LB2R = 1- B2RF ; B1RB2L = 0  } ;

CP25::usage = "Conditional probabilities when r = rSat = 1/n2 = 1/(n1 + 1) = 2/(n1 + n2 + 1), n2 - n1 = 1, k1 and k2 are > zero, w11 = w21 = zero, and w22 and w12 are > zero." ; 
CP25 := { 
     FB1R = (k2 w12 - k1)/(k2 w12 - k1 Q/(1-Q) ) ; 
     B1LF = 1 ; B1RF = FB1R Q/(1-Q) ; B2LF = 0 ; B2RF = 1- FB1R ; 
     FB1L = 1 ; FB2L = 0 ; FB2R = 1- FB1R ; 
     B1LB1R = 0 ; B1RB1L = 0 ; B2LB2R = B1RF ; B2RB2L = FB1R ; 
     B1LB2R = 0 ; B1RB2L = 1 - B1RF  } ;

CP26::usage = "Conditional probabilities when r = rSat = 1/n2 = 1/(n1 + 1) = 2/(n1 + n2 + 1), n2 - n1 = 1, k1 and k2 are > zero, w22 = w21 = zero, and w11 and w12 are > zero." ; 
CP26 := { 
     FB2L = (k2 w12 - k1)/(k2 w12 - k1 Q/(1-Q) ) ; 
     B1LF = 1- FB2L ; B1RF = 0 ; B2LF = FB2L Q/(1-Q) ; B2RF = 1 ; 
     FB1L = 1- FB2L ; FB1R = 0 ; FB2R = 1 ; 
     B1LB1R = FB2L ; B1RB1L = B2LF ; B2LB2R = 0 ; B2RB2L = 0 ; 
     B1LB2R = 0 ; B1RB2L = 1- B2LF  } ;

Calculation of average cluster length

averageClusterLength::usage = "averageClusterLength is the average length (in ligands) of a cluster of bound ligands." ;

averageClusterLength := 
      "infinite" /; Q b2RF + (1 - Q) b1RF == 0; averageClusterLength := 
      (1/(Q b2RF + (1 - Q) b1RF)) /; Q b2RF + (1 - Q) b1RF > 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, printWarning1 ; needPrintWarning1 = False, 
               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} ] ] ; 
     curvePlot = ListPlot [ plotTable, PlotJoined->True, 
          AxesLabel->{" r",yLabel}, PlotLabel->graphLabel, 
          DisplayFunction->Identity, PlotRange->All ] ;
     pointPlot = ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, PlotRange->All ] ; 
     Show [ curvePlot, pointPlot, DisplayFunction->$DisplayFunction, 
          AxesOrigin->{0,yOrigin} ] ; 
     N [ plotTable ] }

makePlot2::usage = "makePlot2 plots the theoretical r/Lf versus r curve (in black) along with corresponding experimental data points (in gray)." ; 
makePlot2 := { 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, 
          printWarning1 ; needPrintWarning1 = False ; 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,1]]}, 
               {i,points2-dropPoints} ] ; 
          AppendTo [ plotTable, 
               {abscissa[[points2+1]], ordinate[[points2+1,1]]} ] , 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,1]]}, {i,points2+1} ] ] ; 
     curvePlot = ListPlot [ plotTable, PlotJoined->True, 
          DisplayFunction->Identity, 
          PlotRange->All ] ; 
     pointPlot = ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, 
          PlotRange->All ] ; 
     plotTable2 = 
          Table [ {rExpArray[[i]], experimentalrLf[[i]]}, {i,points2-1} ]//N; 
     curvePlot2 = ListPlot [ plotTable2, PlotJoined->True, 
          DisplayFunction->Identity, 
          AxesLabel->{"r", "r/Lf"}, 
          PlotLabel->"  Scatchard Plot (with exp. data)", 
          PlotStyle->{GrayLevel [0.75]}, 
          PlotRange->All ] ; 
     pointPlot2 = ListPlot [ plotTable2, 
          PlotStyle->{PointSize[0.015], GrayLevel [0.75] }, 
          DisplayFunction->Identity, 
          PlotRange->All ] ; 
     Show [ curvePlot2, pointPlot2, curvePlot, pointPlot, 
          DisplayFunction->$DisplayFunction, 
          AxesOrigin->{0,0} ] ; 
     N [ plotTable ] }

makeResidualsPlot::usage = "makeResidualsPlot is a plotting routine for the residuals in the plot of the experimental versus theoretical Scatchard plots.  r is plotted on the abscissa." ; 
makeResidualsPlot [ dataList_, curveList_ ] := { 
     residualArray = dataList-curveList ; 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, printWarning2 ] ; 
     xAxisList = Table [ abscissa[[j]], { j, 2, points2} ] ; 
     plotTable = Table [ {xAxisList[[j]], residualArray[[j]]}, 
          {j, points2-dropPoints-1} ] ; 
     curvePlot := ListPlot [ plotTable, PlotJoined->True, 
          AxesLabel->{"r", "residuals"}, 
          PlotLabel->"exp.-theoretical", 
          DisplayFunction->Identity, 
          PlotRange->{{0, rSat}, All} ] ; 
     pointPlot := ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, 
          PlotRange->{{0, rSat}, All} ] ; 
     Show [ curvePlot, pointPlot, 
          DisplayFunction->$DisplayFunction] ; 
     N [ plotTable ] }

makeResidualsPlot2::usage = "makeResidualsPlot2 is a plotting routine for the residuals in the plot of the experimental versus theoretical Scatchard plot.  Point number is plotted on the abscissa." ; 
makeResidualsPlot2 [ dataList_, curveList_ ] := { 
     residualArray = dataList-curveList ; 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, printWarning2 ] ; 
     plotTable = Table [ {j, residualArray[[j]]}, { j, points2-dropPoints-1} ] ; 
     curvePlot := ListPlot [ plotTable, PlotJoined->True, 
          AxesLabel->{"pt. #", "residuals"}, 
          PlotLabel->"exp.-theoretical", 
          DisplayFunction->Identity, 
          PlotRange->{{0, points2-dropPoints-1}, All} ] ; 
     pointPlot := ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, 
          PlotRange->{{0, points2-dropPoints-1}, All} ] ; 
     Show [ curvePlot, pointPlot, 
          DisplayFunction->$DisplayFunction] ; 
     N [ plotTable ] }

makePlot3::usage = "makePlot3 plots the Scatchard plot for asymmetric ligands on an anisotropic lattice (in black) along with the corresponding reference symmetric ligand Scatchard plot (in gray), which is defined only for binding in a single mode (n1 = n2)." ; 
makePlot3 := { 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, printWarning1 ; needPrintWarning1 = False ] ; 
     Print ["w for reference symmetric ligand Scatchard plot = ", w//N ] ; 
     If [ !n==n1 && needReference, 
          Print ["n for reference symmetric ligand Scatchard plot = ", n//N ] ] ; 
     If [ flag, 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,1]]}, 
               {i,points2-dropPoints} ] ; 
          AppendTo [ plotTable, 
               {abscissa[[points2+1]], 0} ] , 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,1]]}, {i,points2+1} ] ] ; 
     curvePlot = ListPlot [ plotTable, PlotJoined->True, 
          AxesLabel->{" r","r/Lf"}, 
          PlotLabel->"Scatchard Plot (with ref. plot)", 
          DisplayFunction->Identity, PlotRange->All ] ; 
     pointPlot = ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, PlotRange->All ] ; 
     If [ needReference,           
          plotTable2 = 
               Table [ {abscissa[[i]], ordinate2[[i]]}, {i,points2+1} ]//N; 
          curvePlot2 = ListPlot [ plotTable2, PlotJoined->True, 
               AxesLabel->{" r","r/Lf"}, 
               PlotLabel->"Scatchard Plot (with ref. plot)", 
               DisplayFunction->Identity, 
               PlotStyle->{GrayLevel [0.75]}, 
               PlotRange->All ] ; 
          pointPlot2 = ListPlot [ plotTable2, 
               PlotStyle->{PointSize[0.015], GrayLevel [0.75] }, 
               DisplayFunction->Identity, PlotRange->All ]  ] ; 
     If [ needReference,           
          Show [ curvePlot2, pointPlot2, curvePlot, pointPlot, 
               DisplayFunction->$DisplayFunction, 
               AxesOrigin->{0,0} ], 
          Show [ curvePlot, pointPlot, 
               DisplayFunction->$DisplayFunction, 
               AxesOrigin->{0,0} ] ] ; 
     N [ plotTable ] }

qPlot::usage = "qPlot is a routine for plotting q1L (the fraction of bound ligands with their heads to the left)." ; 
qPlot := { 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, 
          printWarning2 ; 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,6]]}, 
               {i,points2-dropPoints} ] ; 
          If [ !NEchoice=="NE1" && !NEchoice=="NE2a", 
               AppendTo [ plotTable, 
                    {abscissa[[points2+1]], ordinate[[points2+1,6]]} ] ], 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,6]]}, {i,points2+1} ] ] ; 
     curvePlot = ListPlot [ plotTable, PlotJoined->True, 
          PlotLabel->"ligand orientation", 
          AxesLabel->{" r","   q1L"}, 
          DisplayFunction->Identity, PlotRange->{-0.025,1.025} ] ;
     pointPlot = ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, PlotRange->{-0.025,1.025} ] ; 
     Show [ curvePlot, pointPlot, 
          DisplayFunction->$DisplayFunction ] ; 
     N [ plotTable ] }

latticeCoveragePlot::usage = "latticeCoveragePlot is a routine for plotting the fraction of the lattice covered by ligands." ; 
latticeCoveragePlot:= { 
     latticeCoverageTable = abscissa ; 
     Do [ 
          latticeCoverageTable[[j]] = 
               ordinate[[j,6]] abscissa[[j]] n1 + (1-ordinate[[j,6]]) abscissa[[j]] n2, 
               (* = q r n1 + (1-q) r n2 *) 
          { j, points2+1 } ] ; 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, printWarning2 ; 
          plotTable = Table [ {abscissa[[j]], latticeCoverageTable[[j]]}, 
               {j,points2-dropPoints} ] ; 
          If [ lastPoint, 
               AppendTo [ plotTable, 
                    {abscissa[[points2+1]], latticeCoverageTable[[points2+1]]} ] ], 
          plotTable = Table [ {abscissa[[j]], latticeCoverageTable[[j]]}, { j, points2+1} ] ] ; 
     curvePlot := ListPlot [ plotTable, PlotJoined->True, 
          AxesLabel->{" r", "l.c."}, 
          PlotLabel->"  fraction of lattice covered", 
          DisplayFunction->Identity, 
          PlotRange->{{-0.025 rSat,1.025 rSat},{-0.025,1.025}} ] ; 
     pointPlot := ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, 
          PlotRange->{{-0.025 rSat,1.025 rSat},{-0.025,1.025}} ] ; 
     Show [ curvePlot, pointPlot, DisplayFunction->$DisplayFunction] ; 
     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} ] ] ; 
     curvePlot := ListPlot [ plotTable, PlotJoined->True, 
          AxesLabel->{" r", label}, 
          DisplayFunction->Identity, 
          PlotRange->{{-0.025 rSat,1.025 rSat},{-0.025,1.025}} ] ; 
     pointPlot := ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, 
          PlotRange->{{-0.025 rSat,1.025 rSat},{-0.025,1.025}} ] ; 
     Show [ curvePlot, pointPlot, DisplayFunction->$DisplayFunction] ; 
     N [ plotTable ] }

twoDplot::usage = "twoDplot is a routine for plotting the fractions of clusters of different lengths, or the fractions of bound ligands in clusters of different lengths, at a given value of r." ; 
twoDplot [ yLabel_, maxY_, twoDarray_ ] := { 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, 
          If [ needPrintWarning1, 
               printWarning1 ; needPrintWarning1 = False, 
               printWarning2 ]  ] ; 
     If [ (flag==False) || (rCluster2D rSat < rError) || (rCluster2D==1), 
          twoDTable = Table [ {j, twoDarray[[j]]}, {j, maxLength} ] ; 
          curvePlot = ListPlot [ twoDTable, PlotJoined->True, 
               DisplayFunction->Identity, 
               AxesLabel->{" length", yLabel}, 
               PlotRange->{
                    {1-0.025 (maxLength-1), maxLength+0.025 (maxLength-1)},
                    {-0.025 maxY, 1.025 maxY}} ] ; 
          pointPlot = ListPlot [ twoDTable, PlotStyle->PointSize[0.015], 
               DisplayFunction->Identity, 
               PlotRange->{
                    {1-0.025 (maxLength-1), maxLength+0.025 (maxLength-1)},
                    {-0.025 maxY, 1.025 maxY}} ] ; 
          Show [ curvePlot, pointPlot, 
               DisplayFunction->$DisplayFunction ] ] } ;

rCluster::usage = "rCluster is a number between zero and one, the fractional saturation of the lattice for which the calculations will be performed." ; 
maxLength::usage = "maxLength is the maximum cluster length to be considered." ;

printHeadings::usage = "printHeadings prints the column headings for the cluster length distribution data." ; 
printHeadings := { 
     If [ flag, printWarning2 ; Print [ " " ] ] ; 
     Print[ "r/    Cluster  Fraction (f)  Running     Fraction (f)  Running"] ; 
     Print[ "rSat  Length    Clusters       Sum         Ligands       Sum"] ; 
     Print[ " "] } ;

clusterDistribution::usage = "clusterDistribution is a routine for calculating the fractions of clusters of different lengths, and the fractions of bound ligands in clusters of different lengths, at a given value of r." ; 

(*  This group of nested do-loops calculates the fractions of clusters having different lengths and the fractions of bound ligands residing in clusters of different lengths (using left-to-right conditional probabilities).  arrangement is a number that represents (in binary) the orientations of ligands in a cluster (with zeros and ones representing ligands with their heads to the right and left, respectively).  
     
The probability that a cluster is of a given length is the probability that a free lattice residue is bordered on the right by a cluster of that length divided by the probability that the free lattice residue is bordered on the right by a cluster of any length. The sum of the probabilities for finding clusters of each possible arrangement for a given length, fracClusters, is the probability that a randomly chosen cluster is that long.  
     
The probability that a bound ligand is in a cluster of a given arrangement is the probability that a randomly chosen bound ligand has the orientation appropriate for a given site in the arrangement (Q or 1-Q, depending on the site) times the probabilities of finding all the remaining ligands required by the arrangement (the product is the same no matter what site in the arrangment is initially chosen) times the number of ligands in the cluster.  The sum of the probabilities of finding a ligand in each possible arrangement for a cluster of a given length, fracLigands, is the probability that a randomly chosen bound ligand is in a cluster of that length.  *) 
     
clusterDistribution := { 
(* Find the largest r value <= (rCluster rSat) for which the neighbor-effect values are already calculated *) 
     initial = 1 ; jj =. ; alreadyDone = False ; 
     Do [ r = abscissa[[jj]] ; 
          If [ r < rCluster rSat, initial = jj ] ; 
          If [ r == rCluster rSat, alreadyDone = True ; initial = jj  ] , 
     { jj, 2, points2+1} ] ; 
(* Calculate the conditional probabilities for r = rCluster rSat *) 
     r = rCluster rSat ; 
     If [ r==0, {b1LF,b2LF,b1RF,b2RF,fB1L,fB2L,b1RB1L,b2RB2L,
               b1RB2L,b2RB1L,Q} = 
          {1,1,1,1,0,0,0,0,0,0, k1/(k1 + k2) } ] ; 
     If [ r >0 && r < rSat, 
               E1L = ordinate[[initial,2]] ; E1R = ordinate[[initial,3]] ; 
          E2L = ordinate[[initial,4]] ; E2R = ordinate[[initial,5]] ; 
          Q = ordinate[[initial,6]] ; 
          If [ alreadyDone==False, 
               If [ needCheck, 
                    If [ flag==True && r >= rError, 
                         {E1L,E1R,E2L,E2R,Q} = {1,1,1,1,0.5} ] ; 
                    If [ (flag==False) || (r < rError), ne ] ; 
                    checkLimits, 
                    ne ] ] ; 
          X = r/F ; fF = ff ; 
          {b1LF,b2LF,b1RF,b2RF,fB1L,fB2L,
               b1RB1L,b2RB2L,b1RB2L,b2RB1L} = 
          {fF/E1L, fF/E2L, fF/E1R, fF/E2R, Q X fF/E1L, (1-Q) X fF/E2L,
               Q X w11 fF/(E1R E1L), (1-Q) X w22 fF/(E2R E2L), 
               (1-Q) X w12 fF/(E1R E2L), Q X w21 fF/(E1L E2R)} ] ; 
     If [ r==rSat, 
          {b1LF,b2LF,b1RF,b2RF,fB1L,fB2L,b1RB1L,b2RB2L,
               b1RB2L,b2RB1L,Q} = 
          {B1LF,B2LF,B1RF,B2RF,FB1L,FB2L,B1RB1L,B2RB2L,
               B1RB2L,B1LB2R,ordinate[[points2+1,6]]} ] ; 
(* Calculate the cluster length distribution *) 
     totalClusters = 0 ; maxClusters = 0 ; 
     totalLigands = 0 ; maxLigands = 0 ; 
     clusterLength =. ; 
     Do [ fracClusters = 0 ; fracLigands = 0 ; 
          If [ r==0, 
               If [ clusterLength==1, fracClusters = 1 ; fracLigands = 1 ] ] ; 
          If [ r==rSat && n1==n2, 
               If [ clusterLength==1 && rSat==1/(n1 + 1), 
                    fracClusters = 1 ; fracLigands = 1 ] ; 
               If [ clusterLength==2 && rSat==2/( 2 n1 + 1), 
                    fracClusters = 1 ; fracLigands = 1 ] ] ; 
          If [ r >0 && ( r < rSat || !n1==n2 ) && ((flag==False) || (r < rError)), 
               If [ Q >0, 
                    If [ clusterLength==1, cpProduct = b2RF, 
                         cpProduct = b2RB1L^(clusterLength-1) b2RF ] ; 
                    fracClusters = fB1L cpProduct ; 
                    fracLigands = Q b1LF cpProduct ] ; 
               If [ Q <1, 
                    If [ clusterLength==1, cpProduct = b1RF, 
                         cpProduct = b1RB2L^(clusterLength-1) b1RF ] ; 
                    fracClusters = fracClusters + fB2L cpProduct ; 
                    fracLigands = fracLigands + (1-Q) b2LF cpProduct ] ; 
               If [ Q >0 && Q <1, 
                    arrangement =. ; 
                    Do [ If [ Mod [ arrangement,2]==0, 
                                leftLigand = 1 ; cpProduct = b1RF, 
                                leftLigand = 2 ; cpProduct = b2RF ] ; 
                         arrangement = Quotient [ arrangement, 2 ] ; 
                         cpFactor =. ; 
                         Do [ If [ leftLigand==2, rightLigand = 1, 
                              rightLigand = 2 ] ; 
                              If [ Mod [ arrangement, 2 ]==0, leftLigand = 1, 
                                   leftLigand = 2 ] ; 
                              arrangement = Quotient [ arrangement, 2 ] ; 
                              If [ leftLigand==rightLigand==1, 
                                   cpProduct = cpProduct b1RB1L ] ; 
                              If [ leftLigand==rightLigand==2, 
                                   cpProduct = cpProduct b2RB2L ] ; 
                              If [ leftLigand==1 && rightLigand==2, 
                                   cpProduct = cpProduct b1RB2L ] ; 
                              If [ leftLigand==2 && rightLigand==1, 
                                  cpProduct = cpProduct b2RB1L ] , 
                         { cpFactor, 1, clusterLength-1} ] ; 
                         If [ leftLigand==2, 
                              fracClusters = fracClusters + cpProduct fB1L ; 
                              fracLigands = fracLigands + cpProduct Q b1LF ] ; 
                         If [ leftLigand==1, 
                              fracClusters = fracClusters + cpProduct fB2L ; 
                              fracLigands = fracLigands + cpProduct (1-Q) b2LF ] , 
                    { arrangement, 1, 2^clusterLength-2} ]  ] ; 
               fracClusters = fracClusters/(fB1L + fB2L) ; 
               fracLigands = fracLigands clusterLength ] ; 
          If [ r >0 && r < rSat && flag && r >= rError, 
               fracClusters = I ; fracLigands = I ] ; 
          clusterArray[[clusterLength]] = fracClusters ; 
          totalClusters = totalClusters + fracClusters ; 
          ligandArray[[clusterLength]] = fracLigands ; 
          totalLigands = totalLigands + fracLigands ; 
          If [ ((flag==False) || (r < rError)), 
               If [ fracClusters >maxClusters, maxClusters = fracClusters ] ; 
               If [ fracLigands >maxLigands, maxLigands = fracLigands ] ] ; 
          Print [ rCluster//N, "    ", 
               clusterLength, "      ", 
               NumberForm [ fracClusters//N, 6], "     ", 
               totalClusters//N,"     ", 
               NumberForm [ fracLigands//N, 6 ], "    ", 
               totalLigands//N ], 
     { clusterLength, 1, maxLength} ]  } ;

chooseTicks::usage = "chooseTicks is a routine for selecting the number and position of the tick marks in the lattice saturation axis (r/rSat) in the 3 dimensional plots of cluster length distribution." ; 
chooseTicks := { 
     If [ points3D==1, tickChoice = 1, 
          If [ points3D==2, tickChoice = 2, 
               If [ Mod[points3D,5]==0, tickChoice = 5 ; 
                    tickStep = points3D/5, 
                    If [ Mod[points3D,4]==0, tickChoice = 4 ; 
                         tickStep = points3D/4, 
                         If [ Mod[points3D,3]==0, tickChoice = 3 ; 
                              tickStep = points3D/3, 
                              tickChoice = 6 ] ] ] ] ] ; 
      If [ tickChoice==1, tickArray = {{1,0},{2,1}} ; 
           tickArray0 = {{1,0}} ; tickArray1 = {{2,1}} ] ; 
      If [ tickChoice==2, tickArray = {{1,0},{2,0.5},{3,1}} ; 
           tickArray0 = {{1,0},{2,0.5}} ; tickArray1 = {{2,0.5},{3,1}} ] ; 
      If [ tickChoice==3, 
           tickArray = {{1,0},{1+tickStep,0.33},{1+2 tickStep,0.67}, 
                {1+3 tickStep,1}} ; 
           tickArray0 = {{1,0},{1+tickStep,0.33},{1+2 tickStep,0.67}} ; 
           tickArray1 = {{1+tickStep,0.33},{1+2 tickStep,0.67}, 
                {1+3 tickStep,1}} ] ; 
      If [ tickChoice==4, 
           tickArray = {{1,0},{1+tickStep,0.25},{1+2 tickStep,0.5}, 
                {1+3 tickStep,0.75},{1+4 tickStep,1}} ; 
           tickArray0 = {{1,0},{1+tickStep,0.25},{1+2 tickStep,0.5}, 
                {1+3 tickStep,0.75}} ; 
           tickArray1 = {{1+tickStep,0.25},{1+2 tickStep,0.5}, 
                {1+3 tickStep,0.75},{1+4 tickStep,1}} ] ; 
      If [ tickChoice==5, 
           tickArray = {{1,0},{1+tickStep,0.2},{1+2 tickStep,0.4}, 
                {1+3 tickStep,0.6},{1+4 tickStep,0.8},{1+5 tickStep,1}} ; 
           tickArray0 = {{1,0},{1+tickStep,0.2},{1+2 tickStep,0.4}, 
                {1+3 tickStep,0.6},{1+4 tickStep,0.8}} ; 
           tickArray1 = {{1+tickStep,0.2},{1+2 tickStep,0.4}, 
                {1+3 tickStep,0.6},{1+4 tickStep,0.8},{1+5 tickStep,1}} ] ; 
      If [ tickChoice==6, 
            num = 1 ; 
            While [ num/points3D < 0.1, num = num+1 ] ; 
            num2 = points3D-num ; 
            While [ !Mod[num2,3]==0 && !Mod[num2,4]==0, 
                 num2 = num2-1 ] ; 
            If [ Mod[num2,4]==0, tickChoice2 = 4 ; 
                 tickStep = num2/4, 
                 If [ Mod[num2,3]==0, tickChoice2 = 3 ; 
                      tickStep = num2/3 ] ] ; 
            tickValue = tickStep/points3D ; 
            If [ tickChoice2==3, 
                 tickArray = {{1,0}, 
                      {1+tickStep, N [ Round [ 100 tickValue]/100 ] }, 
                      {1+2 tickStep, N [ Round [ 200 tickValue]/100 ] }, 
                      {1+3 tickStep, N [ Round [ 300 tickValue]/100 ] }, 
                      {points3D+1, 1}} ; 
                 tickArray2 = { 
                      {1+tickStep, N [ Round [ 100 tickValue]/100 ] }, 
                      {1+2 tickStep, N [ Round [ 200 tickValue]/100 ] }, 
                      {1+3 tickStep, N [ Round [ 300 tickValue]/100 ] }} ; 
                 tickArray3 = {{points3D+1-3 tickStep, 
                           N [ Round [ 100 (1-3 tickValue)]/100 ] }, 
                      {points3D+1-2 tickStep, 
                           N [ Round [ 100 (1-2 tickValue)]/100 ] }, 
                      {points3D+1-tickStep, 
                           N [ Round [ 100 (1-tickValue)]/100 ] }} ; 
                 If [ tickStep > (points3D-3 tickStep), 
                      tickArray0 = tickArray2 ; tickArray1 = tickArray3, 
                      tickArray0 = tickArray3 ; tickArray1 = tickArray2 ] ] ; 
            If [ tickChoice2==4, 
                 tickArray = {{1,0}, 
                      {1+tickStep, N [ Round [ 100 tickValue]/100 ] }, 
                      {1+2 tickStep, N [ Round [ 200 tickValue]/100 ] }, 
                      {1+3 tickStep, N [ Round [ 300 tickValue]/100 ] }, 
                      {1+4 tickStep, N [ Round [ 400 tickValue]/100 ] }, 
                      {points3D+1, 1}} ; 
                 tickArray2 = { 
                      {1+tickStep, N [ Round [ 100 tickValue]/100 ] }, 
                      {1+2 tickStep, N [ Round [ 200 tickValue]/100 ] }, 
                      {1+3 tickStep, N [ Round [ 300 tickValue]/100 ] }, 
                      {1+4 tickStep, N [ Round [ 400 tickValue]/100 ] }} ; 
                 tickArray3 = {{points3D+1-4 tickStep, 
                           N [ Round [ 100 (1-4 tickValue)]/100 ] }, 
                      {points3D+1-3 tickStep, 
                           N [ Round [ 100 (1-3 tickValue)]/100 ] }, 
                      {points3D+1-2 tickStep, 
                           N [ Round [ 100 (1-2 tickValue)]/100 ] }, 
                      {points3D+1-tickStep, 
                           N [ Round [ 100 (1-tickValue)]/100 ] }} ; 
                 If [ tickStep > (points3D-4 tickStep), 
                      tickArray0 = tickArray2 ; tickArray1 = tickArray3, 
                      tickArray0 = tickArray3 ; tickArray1 = tickArray2 ] ] ; 
            PrependTo [ tickArray0, {1,0} ] ; 
            AppendTo [ tickArray1, {points3D+1, 1} ] ] } ;

threeDplot::usage = "threeDplot is a routine for plotting the fractions of clusters of different lengths, or the fractions of bound ligands in clusters of different lengths, as a function of r." ; 
threeDplot [ xView_, yView_, zView_, xLabel_, yLabel_, zLabel_, 
     threeDarray_, tickArray_ ] := { 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, 
          If [ needPrintWarning1, 
               printWarning1 ; needPrintWarning1 = False, 
               printWarning2 ] ; Off [ Graphics3D::notnum ] ] ; 
     ListPlot3D [ threeDarray, PlotRange->{0,1}, 
          BoxRatios->{1,1,1}, Lighting -> True, 
          Ticks->{Automatic,tickArray,Automatic}, 
          LightSources -> {{{1,0,1},RGBColor[1,0,0]}, 
               {{1,1,1},RGBColor[0,1,0]}, {{0,1,1},RGBColor[0,0,1]}}, 
          ViewPoint -> {xView, yView, zView}, 
          AxesLabel->{xLabel, yLabel, zLabel} ] ; 
     If [ flag, On [ Graphics3D::notnum ] ]  } ;

Functions for calculation of correlation coefficient (R) for comparison of theoretical and experimental data

\(ssTot[
        list_] := {aveY = Plus@@list\/Length[list]; 
          Plus@@\(\((list - aveY)\)\^2\)}[[1
        ]]; 

    \(ssReg[list1_, list2_] := Plus@@\(\((list1 - list2)\)\^2\); \)

      \(corrCoeff[dataList_, curveList_] := 
        \((1 - ssReg[dataList, curveList]\/ssTot[dataList])\)\^0.5\); \)