The Binding of Asymmetric Ligands to an Infinite Anisotropic Linear Lattice
by Alan R. Wolfe

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,
you should get a dialog box asking whether you want to automatically
evaluate all the initialization cells; click Yes (otherwise, execute the
Action-Evaluate Initialization command).  The need to respond to the dialog
box can be avoided by selecting the "Always" option for "On opening a
Notebook, load initialization cells:" at the bottom of the Action
Preferences box under the Edit menu.  After evaluation, close this notebook
(without quitting Mathematica ) and open the output notebook.  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 ] ;

Establish Mathematica version number

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

Choose whether output is to be saved

saveOutput = False ;

Definitions

points::usage = "points is the number of 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 := (k1+k2) (1-n r) (e^2/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 [ (F/2) 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] >= 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 [ (F/2) (ff^(n1-1) k1 E1L E2R + ff^(n2-1) k2 E1R
E2L) ] ;
             If [ k1 >0 && k2 >0, deltaA = "undefined", deltaA = 0 ] } ;

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-x + Sqrt [(1-x)^2 + 4 x w])/2 ;
     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-X + ( (1-X)^2 +
          4 X (k2 w11 + k1 w21)/(k1 + k2) )^0.5 )/2 } ;

ne3::usage = "Equation for neighbor-effect parameters when r < rSat and
k1 = zero." ;
ne3 := { E1R = E2L = ( 1 - X + ( (1 - X)^2 + 4 X w12 )^0.5 )/2 ; Q = 0 } ;

ne4::usage = "Equation for neighbor-effect parameters when r < rSat and
k2 = zero." ;
ne4 := { E1L = E2R = ( 1 - X + ( (1 - X)^2 + 4 X w21 )^0.5 )/2 ; 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/e1r) (1+ q/(z e1l)+ (1-q)/(z e1r) )^(n2-n1) ==
               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/e1r) (1+ q x/e1l+ (1-q) x/e1r )^(n2-n1) ==
               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/k2) w22)/(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/k1) w22)/(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/k2) w11)/(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/k1) w11)/(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 + (1 + 4 w22 k1/k2)^(-0.5) )/2 ;
     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 + (1 + 4 w22 k2/k1)^(-0.5) )/2 ;
     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 + (1 + 4 w11 k1/k2)^(-0.5) )/2 ;
     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 + (1 + 4 w11 k2/k1)^(-0.5) )/2 ;
     E1R = w11 (2 Q -1)/Q ;
     {E1L,E2L,E2R} = {0,0,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 = (k-w21 + Sqrt [ (k-w21)^2 + 4 w22 k ] )/2 ;
     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 = (k-w21 + Sqrt [ (k-w21)^2 + 4 w11 k ] )/2 ;
     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 = (k-w12 + Sqrt [ (k-w12)^2 + 4 w22 k ] )/2 ;
     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 = (k-w12 + Sqrt [ (k-w12)^2 + 4 w11 k ] )/2 ;
     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 ;

Calculation of the initial slope

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

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,points-dropPoints} ] ;
          If [ lastPoint || j==1 || j==8,
               AppendTo [ plotTable, {abscissa[[points+1]],
                    ordinate[[points+1,j]]} ] ] ,
          plotTable = Table [ {abscissa[[i]], ordinate[[i,j]]},
{i,points+1} ] ] ;
     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,
          If [ versionIsLessThan2, Axes->{0,yOrigin},
               AxesOrigin->{0,yOrigin} ] ] ;
     N [ plotTable ] }

makePlot2::usage = "makePlot2 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)." ;
makePlot2 := {
     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,points-dropPoints} ] ;
          AppendTo [ plotTable,
               {abscissa[[points+1]], 0} ] ,
          plotTable = Table [ {abscissa[[i]], ordinate[[i,1]]},
{i,points+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,points+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,
               If [ versionIsLessThan2, Axes->{0,0}, AxesOrigin->{0,0} ] ],
          Show [ curvePlot, pointPlot,
               DisplayFunction->$DisplayFunction,
               If [ versionIsLessThan2, Axes->{0,0}, 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,points-dropPoints} ] ;
          If [ !NEchoice=="NE1" && !NEchoice=="NE2a",
               AppendTo [ plotTable,
                    {abscissa[[points+1]], ordinate[[points+1,6]]} ] ],
          plotTable = Table [ {abscissa[[i]], ordinate[[i,6]]},
{i,points+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 ] }

makeCPplot::usage = "makeCPplot is a plotting routine for conditional
probabilities." ;
makeCPplot [ cpFunction_, label_ ] := {
     cpTable = Map [ cpFunction , ordinate ] ;
     If [ saveOutput==False,
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ;
     If [ flag, printWarning2 ;
          plotTable = Table [ {abscissa[[j]], cpTable[[j]]},
               {j,points-dropPoints} ] ;
          If [ lastPoint,
               AppendTo [ plotTable,
                    {abscissa[[points+1]], cpTable[[points+1]]} ] ],
          plotTable = Table [ {abscissa[[j]], cpTable[[j]]}, { j, points+1}
] ] ;
     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 ] ] } ;

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.

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 := {
     (* 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, points+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[[points+1,6]]} ] ;
     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 ] ]  } ;