Binding of Symmetric Ligands
to an Infinite Linear Lattice
by Alan R. Wolfe

This is the initialization notebook for the symmetic ligand case.  It needs
to be evaluated only once at the beginning of a Mathematica  session.
After opening this notebook, a dialog box should appear 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. You should 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." ;
w::usage = "w is the cooperativity parameter for interactions between
symmetric ligands bound to a lattice." ;
r::usage = "r is the ratio of bound ligands to lattice residues (plotted
on the abscissa of a Scatchard plot)." ;
x::usage = "x is r/f, where f = 1 - n r is the fraction of lattice residues
that are free." ;
rLf::usage = "rLf is r divided by the free ligand concentration (plotted
on the ordinate of a Scatchard plot)." ;
e::usage = "e is the neighbor-effect parameter for free binding sites." ;
rSat::usage = "rSat is the limiting value of r, i.e., the ratio of ligands
to lattice residues when the lattice is saturated." ;

rLf :=  k (1 - n r)  (1 + x/e )^(1-n) e^2 ;

Expressions for the neighbor-effect parameter when r is below its maximum value

ne::usage = "ne gives the value of the neighbor-effect parameter
when r is below its maximum value." ;
ne :=  If [ w==1, (* THEN *) e = 1,
         (* ELSE *) e := ( 1 - x + ( (1 - x)^2 + 4 x w )^0.5 )/2 ] ;

Determination of the maximum value of r

rSatDetermine::usage = "A routine to determine the ratio of bound
ligands to lattice residues and the values of the conditional
probabilities at at saturation." ;
rSatDetermine := {
     If [ w==0, (* THEN *) rSat = 1/(n + 1); CP2,
          (* ELSE *) rSat = 1/n; CP1 ]  } ;

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." ;
bf::usage = "bf is the probability that a randomly chosen bound
ligand end is bordered by a free lattice residue." ;
fb::usage = "fb is the probability that a randomly chosen free
lattice residue is bordered on a given side by a bound ligand end." ;
bb::usage = "bb is the probability that a randomly chosen bound
ligand end is bordered by another bound ligand end." ;

ff := 1 / (1 + x/e ) ;
bf [h_] := h[[4]]/h[[2]] /; h[[1]] > 0 ; (* = ff/e *)
bf [h_] := bfSat /; h[[1]]==0 ;
fb [h_] := h[[3]] h[[4]]/h[[2]] /; h[[1]] > 0 ; (* = x ff/e *)
fb [h_] := 1 /; h[[1]]==0 ;
bb [h_] := h[[3]] w h[[4]]/h[[2]]^2 /; h[[1]] > 0 ; (* = x w ff/e^2 *)
bb [h_] := bbSAT /; h[[1]]==0 ;

Values of conditional probabilities when r is at its limit

bfSat::usage = "bfSat is the value of bf when the lattice is saturated." ;
bbSAT::usage = "bbSAT is the value of bb when the lattice is saturated." ;

CP1::usage = "CP1 gives conditional probability values when
r =  rSat = 1/n." ;
CP1 := { bfSat = 0 ; bbSAT = 1  } ;

CP2::usage = "CP2 gives conditional probability values when
r =  rSat = 1/(n+1), i.e., when w = zero and the lattice is saturated." ;
CP2 := { bfSat = 1 ; bbSAT = 0  } ;

Calculation of average cluster length

averageClusterLength::usage = "averageClusterLength is the average
length (in ligands) of a cluster of bound ligands." ;
averageClusterLength := "infinite" /;  bfCluster == 0 ;
averageClusterLength := 1/bfCluster /;  bfCluster > 0 ;

Calculation of the initial slope

initialSlope := k (2 w - 2 n - 1) ;

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

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 ] ] ;
     plotTable = Table [ {abscissa[[j]], cpTable[[j]]}, { j, points+1} ] ;
     curvePlot := ListPlot [ plotTable, PlotJoined->True,
          AxesLabel->{" r", label},
          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 ] }

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." ;

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

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

clusterDistribution := {
     (* Calculate the conditional probabilities for
         r = rCluster rSat *)
     r = rCluster rSat ;
     If [ r < rSat, x = r/(1 - n r) ;
               {bfCluster, bbCLUSTER} = {ff/e, x w ff/e^2},
          {bfCluster, bbCLUSTER} = {bfSat,bbSAT} ] ;
     totalClusters = 0 ; maxClusters = 0 ;
     totalLigands = 0 ; maxLigands = 0 ;
     clusterLength =. ;
     Map [ Function [ clusterLength,
          If [  bbCLUSTER==0 && clusterLength==1, fracClusters = 1,
                 fracClusters =  bbCLUSTER^(clusterLength-1) bfCluster ] ;
          clusterArray[[clusterLength]] = fracClusters ;
          If [ fracClusters >maxClusters, maxClusters = fracClusters ] ;
          totalClusters = totalClusters + fracClusters ;
          fracLigands = clusterLength fracClusters bfCluster ;
          ligandArray[[clusterLength]] = fracLigands ;
          If [ fracLigands >maxLigands, maxLigands = fracLigands ] ;
          totalLigands = totalLigands + fracLigands ;
          Print [ rCluster//N, "    ",
               clusterLength, "      ",
               NumberForm [ fracClusters//N, 6], "     ",
               totalClusters//N,"     ",
               NumberForm [ fracLigands//N, 6 ], "    ",
               totalLigands//N ] ],
     clusterLengthArray ] } ;

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