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

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

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." ;
k::usage = "k is the binding constant for ligand-lattice interactions." ;
n::usage = "n is the number of lattice residues occupied per bound
ligand." ;
nRef::usage = "nRef is the effective value of n for the reference
symmetric ligand Scatchard plot (nRef = n except when w3 and one
of w1 or w2 are zero, in which case nRef = n + 1/2)." ;
w::usage = "w is the value of the cooperativity parameter for the
reference symmetric ligand Scatchard plot." ;
w1::usage = "w1 is the cooperativity parameter for head-head (b1-b1)
interactions between ligands bound to an isotropic lattice." ;
w2::usage = "w2 is the cooperativity parameter for tail-tail (b2-b2)
interactions between ligands bound to an isotropic lattice." ;
w3::usage = "w3 is the cooperativity parameter for head-tail (b1-b2)
or tail-head interactions between ligands bound to an isotropic lattice." ;
r::usage = "r is the ratio of bound ligands to lattice residues (plotted
on the abscissa of a Scatchard plot)." ;
f::usage = "f = 1 - n r is the fraction of lattice residues that are free." ;
x::usage = "x is r/f." ;
y::usage = "y is x/2, or r/(2 f)." ;
rLf::usage = "rLf is r divided by the free ligand concentration (plotted
on the ordinate of a Scatchard plot)." ;
E1::usage = "E1 is the neighbor-effect parameter for the head (b1) end
of a binding site on an isotropic lattice." ;
E2::usage = "E2 is the neighbor-effect parameter for the tail (b2) end
of a binding site on an isotropic lattice." ;
eref::usage = "eref is the neighbor-effect parameter for the reference
symmetric ligand Scatchard plot." ;
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 - | w1 - w2 |.
(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." ;

rLf := Re [ k f  ff^(n-1) E1 E2 ] ;
rLfRef := k (1 - nRef r) eref^2 (1 + xREF/eref)^(1-nRef) ;

Determination of the maximum value of r

rSatDetermine::usage = "A routine to determine the ratio of bound
ligands to lattice residues at saturation." ;
rSatDetermine := {
     If [ w1==w2==w3==0, rSat = 1/(n + 1),
          If [ w3 == 0 && Xor [ w1==0, w2 == 0 ], rSat = 2/(2 n + 1),
               rSat = 1/n ] ] } ;

Equations for calculating values of neighbor-effect parameters; parameter
values for special cases

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 = (w1 + w2 + 2 w3)/4 ;
     If [ Abs [deltaA] <10^(-12) || !referencePlot,
          needReference = False,
          needReference = True ; nRef = n ;
               If [ w3==0 && Xor [ w1==0, w2==0 ],
                    w = w + 0.5 ; nRef = nRef + 0.5 ] ] ;
     If [ w1==w2==w3==1, ne := ne3 ; neChoice = "ne3",
          If [ w1==w2, ne := ne2 ; neChoice = "ne2",
               ne := ne1 ; neChoice = "ne1" ] ] ;
     If [ neChoice=="ne1", needCheck = True, needCheck = False ] } ;

setLimits::usage = "setLimits is a routine that determines the
largest and smallest physically possible values for the neighbor-effect
parameters and conditional probabilities." ;
setLimits := {
     neTol = 0.0000001 ;
     E1min = Min [ 1, w1, w3 ] - neTol ;
     E1max = Max [ 1, w1, w3 ] + neTol ;
     E2min = Min [ 1, w2, w3 ] - neTol ;
     E2max = Max [ 1, w2, w3 ] + neTol ;
     cpTol = 0.000001 ;
     cpMin = - cpTol ;
     cpMax = 1 + cpTol } ;

checkLimits::usage = "checkLimits is a routine that checks whether
the calculated values of the neighbor-effect parameters and conditional
probabilities are within their permitted ranges, and sets flag = True
(indicating that a computational error has occurred) if one of them is
outside its permitted range." ;
checkLimits := {
     If [ (flag==False) || (r < rError),
          If [ Head [ E1 ] == Complex || Head [ E2 ] == Complex,
               {E1,E2} = Map [ Re, {E1,E2} ] ;
               rError = r ;
               If [ flag==False, needPrintWarning1 = True ] ;
               flag = True ]  ] ;
     If [ (flag==False) || (r < rError),
          If [ E1<E1min || E1>E1max || E2<E2min || E2>E2max,
               rError = r ;
               If [ flag==False, needPrintWarning1 = True ] ;
               flag = True ]  ] ;
     If [ (flag==False) || (r < rError),
          fF = ff ; b1fTest = fF/E1 ; b2fTest = fF/E2 ;
          fb1Test = y fF/E1 ; fb2Test = y fF/E2 ;
          b1b1Test = y w1 fF/E1^2 ; b2b2Test = y w2 fF/E2^2 ;
          b1b2Test = y w3 fF/(E1 E2) ;
          If [ fF < cpMin || fF > cpMax ||
               b1fTest < cpMin || b1fTest > cpMax ||
               b2fTest < cpMin || b2fTest > cpMax ||
               fb1Test < cpMin || fb1Test > cpMax ||
               fb2Test < cpMin || fb2Test > cpMax ||
               b1b1Test < cpMin || b1b1Test > cpMax ||
               b2b2Test < cpMin || b2b2Test > cpMax ||
               b1b2Test < cpMin || b1b2Test > 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 = "ne1 gives equations for the neighbor-effect parameters
when r is below its maximum value (used when w1 does not = w2)." ;
ne1 := {
     {E1,E2} = {e1,e2}/.FindRoot [ {
          e1 (1/y + 1/e1 + 1/e2) == 1/y + w1/e1 + w3/e2 ,
          e2 (1/y + 1/e1 + 1/e2) == 1/y + w3/e1 + w2/e2  } ,
          {e1,E1}, {e2,E2}, AccuracyGoal ->7,
          MaxIterations ->200, WorkingPrecision ->18 ]  } ;

ne1::usage = "ne1 gives equations for the neighbor-effect parameters
when r is below its maximum value (used when w1 does not = w2)." ;
(* Alternate Version *)
ne1 := {
     {E1,E2} = {e1,e2}/.FindRoot [ {
          e1 == (1/y + w1/e1 + w3/e2)/(1/y + 1/e1 + 1/e2) ,
          e2 == (1/y + w3/e1 + w2/e2)/(1/y + 1/e1 + 1/e2)  } ,
          {e1,E1}, {e2,E2}, AccuracyGoal ->7,
          MaxIterations ->200, WorkingPrecision ->18 ]  } ;

ne2::usage = "ne2 gives an equation for the neighbor-effect
parameters when r is below its maximum value and w1 = w2
(used when w1 and w3 are not both 1)." ;
ne2 := E1 = E2 = ( 1 - x + ( (1 - x)^2 + 2 x (w1 + w3) )^0.5 )/2 ;

ne3::usage = "ne3 gives values of 1 for the neighbor-effect parameters
when r is below its maximum value and w1 = w2 = w3 = 1." ;
ne3 := E1 = E2 = 1 ;

startPrint::usage = "startPrint starts the printout of the Scatchard plot." ;
startPrint := {
     Print [ "r           r/Lf" ] ;
     Print [ "0.0         ", k ] } ;

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 - nRef r) ;
               eref = (1-xREF + Sqrt [(1-xREF)^2 + 4 xREF w] )/2 ;
               rLfRef ], abscissa ] ;
          AppendTo [ ordinate2, 0 ]  ] ;
     AppendTo [ abscissa, rSat ] ;
     Print [ rSat//N, "     0" ]  } ;

neSat::usage = "neSat gives equations for neighbor-effect parameters
when r is at its maximum value." ;
neSat := {
     If [ w1==w2 , (* THEN *)
          E1 = E2 = (w1+w3)/2 , (* ELSE *)
          E1 = ( w1 (w3-w2) + (w1-w3) (w1 w2)^0.5 ) / (w1-w2) ;
          E2 = ( w2 (w3-w1) + (w2-w3) (w1 w2)^0.5 ) / (w2-w1) ]  } ;

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." ;
b1f::usage = "b1f is the probability that a randomly chosen
b1 end of a bound ligand is bordered by a free lattice residue." ;
b2f::usage = "b2f is the probability that a randomly chosen
b2 end of a bound ligand is bordered by a free lattice residue." ;
fb1::usage = "fb1 is the probability that a randomly chosen free
lattice residue is bordered on a given side by a b1 end of a bound ligand." ;
fb2::usage = "fb2 is the probability that a randomly chosen free
lattice residue is bordered on a given side by a b2 end of a bound ligand." ;
b1b1::usage = "b1b1 is the probability that a randomly chosen b1 end
of a bound ligand is bordered by another b1 end of a bound ligand." ;
b2b2::usage = "b2b2 is the probability that a randomly chosen b2 end
of a bound ligand is bordered by another b2 end of a bound ligand." ;
b1b2::usage = "b1b2 (= b2b1) is the probability that a randomly chosen
b1 end of a bound ligand is bordered by a b2 end of a bound ligand." ;

ff := 1 / (1 + y/E1 + y/E2) ;
b1f [ h_] :=  h[[5]]/ h[[2]] /;  h[[1]] > 0 ; (* = ff/E1 *)
b1f [ h_] := b1fSat /;  h[[1]]==0 ;
b2f [ h_] :=  h[[5]]/ h[[3]] /;  h[[1]] > 0 ; (* = ff/E2 *)
b2f [ h_] := b2fSat /;  h[[1]]==0 ;
fb1 [ h_] :=  h[[4]]  h[[5]]/ h[[2]] /;  h[[1]] > 0 ; (* = y ff/E1 *)
fb1 [ h_] := fb1Sat /;  h[[1]]==0 ;
fb2 [ h_] :=  h[[4]]  h[[5]]/ h[[3]] /;  h[[1]] > 0 ; (* = y ff/E2 *)
fb2 [ h_] := fb2Sat /;  h[[1]]==0 ;
b1b1 [ h_] :=  h[[4]] w1  h[[5]]/ h[[2]]^2 /;  h[[1]] > 0 ; (* = y w1
ff/E1^2 *)
b1b1 [ h_] := b1b1Sat /;  h[[1]]==0 ;
b2b2 [ h_] :=  h[[4]] w2  h[[5]]/ h[[3]]^2 /;  h[[1]] > 0 ; (* = y w2
ff/E2^2 *)
b2b2 [ h_] := b2b2Sat /;  h[[1]]==0 ;
b1b2 [ h_] :=  h[[4]] w3  h[[5]]/( h[[2]]  h[[3]]) /;  h[[1]] > 0 ;
     (* = y w3 ff/(E1 E2) *)
b1b2 [ h_] := b1b2Sat /;  h[[1]]==0 ;

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

b1fSat::usage = "b1fSat is the value of b1f when the lattice is saturated." ;
b2fSat::usage = "b2fSat is the value of b2f when the lattice is saturated." ;
fb1Sat::usage = "fb1Sat is the value of fb1 when the lattice is saturated." ;
fb2Sat::usage = "fb2Sat is the value of fb2 when the lattice is saturated." ;
b1b1Sat::usage = "b1b1Sat is the value of b1b1 when the lattice is
saturated." ;
b2b2Sat::usage = "b2b2Sat is the value of b2b2 when the lattice is
saturated." ;
b1b2Sat::usage = "b1b2Sat is the value of b1b2 when the lattice is
saturated." ;

cpSat::usage = "A routine to determine the values of the conditional
probabilities when r =  rSat." ;
cpSat := {
     If [ rSat == 1/n && !( Xor [ w1==0,w2==0 ] ), CP1 ;
          cpChoice = "CP1" ] ;
     If [ rSat == 1/n && w1==0 && w2 >0, CP2 ;
          cpChoice = "CP2" ] ;
     If [ rSat == 1/n && w1 >0 && w2==0, CP3 ;
          cpChoice = "CP3" ] ;
     If [ rSat == 1/(n+1), CP4 ; cpChoice = "CP4" ] ;
     If [ rSat == 2/(2n+1),
          If [ w2==0, CP5 ; cpChoice = "CP5",
          CP6 ; cpChoice = "CP6" ] ] } ;

CP1::usage = "CP1 gives expressions for conditional probabilities when
r =  rSat = 1/n and w1 and w2 are both > zero or are both = zero." ;
CP1 := { b1fSat = 0 ; b2fSat = 0 ;
            fb1Sat = 1/((1/E1+ 1/E2 ) E1) ;
            fb2Sat = 1/((1/E1+ 1/E2 ) E2) ;
            b1b1Sat = w1/((1/E1+ 1/E2 ) E1^2) ;
            b2b2Sat = w2/((1/E1+ 1/E2 ) E2^2) ;
            b1b2Sat = w3/((1/E1+ 1/E2 ) E1 E2)  } ;

CP2::usage = "CP2 gives values of conditional probabilities when
r =  rSat = 1/n, w1 = zero, and w2 > zero." ;
CP2 := { b1fSat = 0 ; b2fSat = 0 ; fb1Sat = 1 ; fb2Sat = 0 ;
            b1b1Sat = 0 ; b2b2Sat = 0 ; b1b2Sat = 1  } ;

CP3::usage = "CP3 gives values of conditional probabilities when
r =  rSat = 1/n, w1 > zero, and w2 = zero." ;
CP3 := { b1fSat = 0 ; b2fSat = 0 ; fb1Sat = 0 ; fb2Sat = 1 ;
            b1b1Sat = 0 ; b2b2Sat = 0 ; b1b2Sat = 1  } ;

CP4::usage = "CP4 gives values of conditional probabilities when
r =  rSat = 1/(n+1), i.e. when w1 = w2 = w3 = zero." ;
CP4 := { b1fSat = 1 ; b2fSat = 1 ; fb1Sat = 0.5 ; fb2Sat = 0.5 ;
           b1b1Sat = 0 ; b2b2Sat = 0 ; b1b2Sat = 0  } ;

CP5::usage = "CP5 gives values of conditional probabilities when
r =  rSat = 2/(2n+1) and w2 = w3 = zero." ;
CP5 := { b1fSat = 0 ; b2fSat = 1 ; fb1Sat = 0 ; fb2Sat = 1 ;
            b1b1Sat = 1 ; b2b2Sat = 0 ; b1b2Sat = 0  } ;

CP6::usage = "CP6 gives values of conditional probabilities when
r =  rSat = 2/(2n+1) and w1 = w3 = zero." ;
CP6 := { b1fSat = 1 ; b2fSat = 0 ; fb1Sat = 1 ; fb2Sat = 0 ;
            b1b1Sat = 0 ; b2b2Sat = 1 ; b1b2Sat = 0  } ;

Calculation of average cluster length

averageClusterLength::usage = "averageClusterLength is the average
length (in ligands) of a cluster of bound ligands." ;
averageClusterLength := "infinite" /;
     b1fCluster + b2fCluster == 0 ;
averageClusterLength := 2/(b1fCluster + b2fCluster) /;
     b1fCluster + b2fCluster > 0 ;

Calculation of the initial slope

initialSlope := k (w1/2 + w2/2 + w3 - 2 n - 1) ;

Calculation of the dimensionless quantity deltaA

deltaA := - Abs [ w1-w2 ] ;

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} ] ;
          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 isotropic lattice (in black) along with the
corresponding reference symmetric ligand Scatchard plot (in gray)." ;
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 [ !nRef==n && needReference,
          Print ["n for reference symmetric ligand Scatchard plot = ",
               nRef//N ]] ;
     If [ flag,
          plotTable = Table [ {abscissa[[i]], ordinate[[i,1]]},
               {i,points-dropPoints} ] ;
          AppendTo [ plotTable,
               {abscissa[[points+1]], ordinate[[points+1,1]]} ] ,
          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 ] }

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} ] ;
          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 (always 1/2 on
an isotropic lattice) 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, {b1fCluster,b2fCluster,fb1Cluster,fb2Cluster,
          b1b1Cluster,b2b2Cluster,b1b2Cluster} = {1,1,0,0,0,0,0} ] ;
     If [ r >0 && r < rSat,
          E1 = ordinate[[initial,2]] ; E2 = ordinate[[initial,3]] ;
          f = 1 - n r ; x = r/f ; y = x/2 ;
          If [ alreadyDone==False,
               If [ needCheck,
                    If [ flag==True && r >= rError, {E1,E2} = {1,1} ] ;
                    If [ (flag==False) || (r < rError), ne ] ;
                    checkLimits,
                    ne ] ] ;
          fF = ff ;
          {b1fCluster,b2fCluster,fb1Cluster,fb2Cluster,
               b1b1Cluster,b2b2Cluster,b1b2Cluster} =
               {fF/E1, fF/E2, y fF/E1, y fF/E2, y w1 fF/E1^2,
                    y w2 fF/E2^2, y w3 fF/(E1 E2) } ] ;
     If [ r==rSat, {b1fCluster,b2fCluster,fb1Cluster,fb2Cluster,
          b1b1Cluster,b2b2Cluster,b1b2Cluster} =
          {b1fSat,b2fSat,fb1Sat,fb2Sat,b1b1Sat,b2b2Sat,b1b2Sat} ] ;
     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,
               If [ clusterLength==1 && rSat==1/(n + 1),
                    fracClusters = 1 ; fracLigands = 1 ] ;
               If [ clusterLength==2 && rSat==2/( 2 n + 1),
                    fracClusters = 1 ; fracLigands = 1 ] ] ;
          If [ r >0 && r < rSat && ((flag==False) || (r < rError)),
               If [ clusterLength==1, cpProduct = 2 b2fCluster,
                    cpProduct = 2 b1b2Cluster^(clusterLength-1) b2fCluster ] ;
               fracClusters = cpProduct fb1Cluster ;
               fracLigands = cpProduct b1fCluster ;
               arrangement =. ;
               Do [ If [ Mod [ arrangement, 2]==0,
                           leftLigand = 1 ; cpProduct = b1fCluster,
                           leftLigand = 2 ; cpProduct = b2fCluster ] ;
                    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 b1b1Cluster ] ;
                         If [ leftLigand==rightLigand==2,
                              cpProduct = cpProduct b2b2Cluster ] ;
                         If [ !leftLigand==rightLigand,
                              cpProduct = cpProduct b1b2Cluster ] ,
                    { cpFactor, 1, clusterLength-1} ] ;
                    If [ leftLigand==2,
                         fracClusters = fracClusters + cpProduct fb1Cluster ;
                         fracLigands = fracLigands + cpProduct b1fCluster ] ;
                    If [ leftLigand==1,
                         fracClusters = fracClusters + cpProduct fb2Cluster ;
                         fracLigands = fracLigands + cpProduct b2fCluster ] ,
               { arrangement, 1, 2^clusterLength-2} ] ;
               fracClusters = fracClusters/(fb1Cluster + fb2Cluster) ;
               fracLigands = fracLigands clusterLength/2 ] ;
          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 ] ]  } ;