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

This is the initialization notebook for the isotropic lattice, asymmetic ligand case.  It needs to be evaluated only once per Mathematica  session.  After opening this notebook, execute the "Kernel-> Evaluation-> Evaluate Initialization" command.  Then close this notebook (without quitting Mathematica ) and open the corresponding output notebook.  It is best to  quit Mathematica  and start a new session when switching to a different case.

Turn off General::spell1 warning

Off [ General::spell1 ] ;

Choose whether output is to be saved

saveOutput = False ;

Definitions

points::usage = "points is the number of 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));

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];

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/2) (1 - x + ((((((1 - x))^2) + 2 x (w1 + w3)))^0.5))
              ;

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, E1 = E2 = ((w1 + w3)/2), 
          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;

Define plotting routines

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

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

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

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

makePlot3::usage = "makePlot3 plots the Scatchard plot for asymmetric ligands on an isotropic lattice (in black) along with the corresponding reference symmetric ligand Scatchard plot (in gray)." ; 
makePlot3 := { 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, 
          If [ needPrintWarning1, 
               printWarning1 ; needPrintWarning1 = False, 
               printWarning2 ] ] ; 
     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,points2-dropPoints} ] ; 
          AppendTo [ plotTable, 
               {abscissa[[points2+1]], ordinate[[points2+1,1]]} ] , 
          plotTable = Table [ {abscissa[[i]], ordinate[[i,1]]}, {i,points2+1} ] ] ; 
     curvePlot = ListPlot [ plotTable, PlotJoined->True, 
          AxesLabel->{" r","r/Lf"}, 
          PlotLabel->"Scatchard Plot (with ref. plot)", 
          DisplayFunction->Identity, PlotRange->All ] ; 
     pointPlot = ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, PlotRange->All ] ; 
     If [ needReference, 
          plotTable2 = 
               Table [ {abscissa[[i]], ordinate2[[i]]}, {i,points2+1} ]//N ; 
          curvePlot2 = ListPlot [ plotTable2, PlotJoined->True, 
               AxesLabel->{" r","r/Lf"}, 
               PlotLabel->"Scatchard Plot (with ref. plot)", 
               DisplayFunction->Identity, 
               PlotStyle->{GrayLevel [0.75]}, 
               PlotRange->All ] ; 
          pointPlot2 = ListPlot [ plotTable2, 
               PlotStyle->{PointSize[0.015], GrayLevel [0.75] }, 
               DisplayFunction->Identity, PlotRange->All ]  ] ; 
     If [ needReference,                
          Show [ curvePlot2, pointPlot2, curvePlot, pointPlot, 
               DisplayFunction->$DisplayFunction, 
               AxesOrigin->{0,0} ], 
          Show [ curvePlot, pointPlot, 
               DisplayFunction->$DisplayFunction, 
               AxesOrigin->{0,0} ] ] ; 
     N [ plotTable ] }

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

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

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

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

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

(*  This group of nested do-loops calculates the fractions of clusters having different lengths and the fractions of bound ligands residing in clusters of different lengths (using left-to-right conditional probabilities).  arrangement is a number that represents (in binary) the orientations of ligands in a cluster (with zeros and ones representing ligands with their heads to the right and left, respectively).  
     
The probability that a cluster is of a given length is the 
 probability that a free lattice residue is bordered on the right by a cluster of that length divided by the probability that the free lattice residue is bordered on the right by a cluster of any length.  The sum of the probabilities for finding clusters of each possible arrangement for a given length, fracClusters, is the probability that a randomly chosen cluster is that long.  
     
The probability that a bound ligand is in a cluster of a given arrangement is the probability that a randomly chosen bound ligand has the orientation appropriate for a given site in the arrangement (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.  *)
     
clusterDistribution := { 
(* Find the largest r value <= (rCluster rSat) for which the neighbor-effect values are already calculated *) 
     initial = 1 ; jj =. ; alreadyDone = False ; 
     Do [ r = abscissa[[jj]] ; 
          If [ r < rCluster rSat, initial = jj ] ; 
          If [ r == rCluster rSat, alreadyDone = True ; initial = jj  ] , 
     { jj, 2, points2+1} ] ; 
(* Calculate the conditional probabilities for r = rCluster rSat *) 
     r = rCluster rSat ; 
     If [ r==0, {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} ] ; 
(* Calculate the cluster length distribution *) 
     totalClusters = 0 ; maxClusters = 0 ; 
     totalLigands = 0 ; maxLigands = 0 ; 
     clusterLength =. ; 
     Do [ fracClusters = 0 ; fracLigands = 0 ; 
          If [ r==0, 
               If [ clusterLength==1, fracClusters = 1 ; fracLigands = 1 ] ] ; 
          If [ r==rSat, 
               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 ] ]  } ;

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

ssTot[
        list_] := {aveY = Plus@@(list/Length)[list]; 
          Plus@@(((list - aveY))^2)}[[1
        ]]; ssReg[list1_, list2_] := Plus@@(((list1 - list2))^2); corrCoeff[dataList_, curveList_] := 
      (((1 - ssReg[dataList, curveList]/ssTot[dataList]))^0.5);