Binding of Symmetric Ligands 
to an Infinite 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 symmetic ligand case.  It needs to be evaluated only once at the beginning of a 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.  You should 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." ; 
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, e = 1, 
        e := (1/2) (1 - x + ((((((1 - x))^2) + 4 x w))^0.5))];

Calculation of the initial slope

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

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, rSat = (1/(n + 1)); CP2, 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 ;

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,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 ] ] ; 
     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", 
          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 ] ] ; 
     xAxisList = Table [ abscissa[[j]], { j, 2, points2} ] ; 
     plotTable = Table [ {xAxisList[[j]], residualArray[[j]]}, 
          {j, points2-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 experimental versus theoretical Scatchard plots.  Point number is plotted on the abscissa." ; 
makeResidualsPlot2 [ dataList_, curveList_ ] := { 
     residualArray = dataList-curveList ; 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     plotTable = Table [ {j, residualArray[[j]]}, { j, points2-1} ] ; 
     curvePlot := ListPlot [ plotTable, PlotJoined->True, 
          AxesLabel->{"pt. #", "residuals"}, 
          PlotLabel->"exp.-theoretical", 
          DisplayFunction->Identity, 
          PlotRange->{{0, points2-1}, All} ] ; 
     pointPlot := ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
          DisplayFunction->Identity, 
          PlotRange->{{0, points2-1}, All} ] ; 
     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 ] ] ; 
     plotTable = Table [ {abscissa[[j]], cpTable[[j]]}, { j, points2+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::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." ; 
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} ] ; 
     (* Calculate the cluster length distribution *) 
     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} ]  } ;

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