The Binding of Symmetric Ligands of
Two Kinds to an Infinite Isotropic Linear Lattice - 
Comparison of Experimental and Theoretical Data  
by Alan R. Wolfe
Department of Biopharmaceutical Sciences
University of California
San Francisco, CA  94143-0446

This is the initialization notebook for the case of competitive binding of 2 kinds of symmetic ligand to an isotropic lattice.  This version was expressly designed to permit comparison of experimental and theoretical data.  The experimental data must be in the form of r2 as a function of either total or free ligand 1 concentration.  The notebook 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.  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 positive ligand 1 concentration values for which parameter values will be calculated.  It also equals the number of intervals on the x-axis of the plots." ; 
maxConcL1::usage = "maxConcL1 is the maximum concentration of ligand 1 for which parameter values will be calculated.  This concentration is expressed in the same units as kL1.  This is a concentration of total ligand 1 for type 1 isotherms and a concentration of free ligand 1 for type 2 isotherms." ; 
kL1::usage = "kL1 is the binding constant for ligand 1-lattice interactions." ; 
kL2::usage = "kL2 is the binding constant for ligand 2-lattice interactions." ; 
nL1::usage = "nL1 is the number of lattice residues occupied per bound ligand 1." ; 
nL2::usage = "nL2 is the number of lattice residues occupied per bound ligand 2." ; 
wL1::usage = "wL1 is the cooperativity parameter for ligand 1-ligand 1 interactions on an isotropic lattice." ; 
wL2::usage = "wL2 is the cooperativity parameter for ligand 2-ligand 2 interactions on an isotropic lattice." ; 
wL12::usage = "wL12 is the cooperativity parameter for ligand 1-ligand 2 (or vice versa) interactions on an isotropic lattice." ; 
r1::usage = "r1 is the ratio of bound ligand 1 to lattice residues (plotted on the abscissa of a Scatchard plot).  This variable is stored as the last value in the sets of values in the array named ordinate." ; 
R1::usage = "R1 is the ratio of bound ligand 1 to lattice residues (used only in equations in which it is being solved for)." ; 
abscissa::usage = "The values of the concentration of ligand 1 (total concentration for type=1, free concentration for type=2) are stored in the array named abscissa." ; 
r2::usage = "r2 is the ratio of bound ligand 2 to lattice residues (plotted on the abscissa of a Scatchard plot)." ; 
R2::usage = "R2 is the ratio of bound ligand 2 to lattice residues (used only in equations in which it is being solved for)." ; 
f::usage = "f is the fraction of lattice residues that are free." ; 
F::usage = "F is the fraction of lattice residues that are free (used in equations being solved for R1 and R2)." ; 
x1::usage = "x1 is r1/f." ; 
x2::usage = "x2 is r2/f." ; 
r1L1f::usage = "r1L1f is r1 divided by the free concentration of ligand 1 (plotted on the ordinate of a Scatchard plot)." ; 
r2L2f::usage = "r2L2f is r2 divided by the free concentration of ligand 2 (plotted on the ordinate of a Scatchard plot)." ; 
EL1::usage = "EL1 is the neighbor-effect parameter for a ligand 1 (symmetric) binding site on an isotropic lattice." ; 
EL2::usage = "EL2 is the neighbor-effect parameter for a ligand 2 (symmetric) binding site on an isotropic lattice." ; 
r1Sat::usage = "r1Sat is the limiting value of r1 for a given value of r2, i.e., the ratio of ligand 1 to lattice residues when the lattice is saturated." ; 
r2Sat::usage = "r2Sat is the limiting value of r2, i.e., the ratio of ligand 2 to lattice residues when the lattice is saturated (and no other ligand is present)." ; 
el2::usage = "el2 is the neighbor-effect parameter for ligand 2 (used only in calculating r2 on a lattice without any bound ligand 1)." ; 
fracRSat2::usage = "fracRSat2 is the fraction of the maximum possible ratio of bound ligand 2 to lattice residues that is present on the lattice." ; 
freeL1::usage = "freeL1 is the concentration of free ligand 1." ; 
totalL1::usage = "totalL1 is the total concentration of ligand 1." ; 
freeL2::usage = "freeL2 is the concentration of free ligand 2." ; 
totalL2::usage = "totalL2 is the total concentration of ligand 2." ; 
polymer::usage = "polymer is the total lattice residue concentration." ;

r1L1f := Re[kL1 f (ff^(nL1 - 1)) (EL1^2)]; r2L2f := Re[kL2 f (ff^(nL2 - 1)) (EL2^2)]; latticeCoverage := r1 nL1 + r2 nL2;

Equations for calculating values of neighbor-effect parameters on an unsaturated lattice

neChoose::usage = "A routine to determine which set of equations for neighbor-effect parameters should be used when r1 < r1Sat." ; 
neChoose := { 
     If [ wL1==wL2==wL12, 
          If [ wL1==wL2==wL12==1, 
               ne := ne5 ; neChoice = "ne5", ne := ne4 ; neChoice = "ne4" ], 
          ne := ne3 ; neChoice = "ne3" ] ; 
     needCheck = True } ;

setLimits::usage = "setLimits is a routine that determines the largest and smallest physically possible values for the neighbor-effect parameters." ; 
setLimits := { 
     neTol = 0.0000001 ; 
     r2min =  - neTol ; 
     r2max = r2Sat + neTol ; 
     EL1min = Min [ 1, wL1, wL12 ] - neTol ; 
     EL1max = Max [ 1, wL1, wL12 ] + neTol ; 
     EL2min = Min [ 1, wL2, wL12 ] - neTol ; 
     EL2max = Max [ 1, wL2, wL12 ] + neTol ; 
     cpTol = 0.000001 ; 
     cpMin = - cpTol ; 
     cpMax = 1 + cpTol ; 
     r1min =  - neTol ; 
     r1max = r1Sat + neTol } ;

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, 
          If [ Head [ EL1] == Complex || Head [ EL2] == Complex || 
               Head [ r2] == Complex, 
               {EL1,EL2,r2} = Map [ Re, {EL1,EL2,r2} ] ; 
               r1Error = r1 ; flag = True ]  ] ; 
     If [ flag == False, 
          If [ r1<r1min || r1>r1max || r2<r2min || r2>r2max || 
               EL1<EL1min || EL1>EL1max || 
               EL2<EL2min || EL2>EL2max || f<0 || f>1, 
               r1Error = r1 ; flag = True ]  ] ; 
     If [ flag==False, 
          fF = ff ; bL1fTest = fF/EL1 ; bL2fTest = fF/EL2 ; 
          fbL1Test = x1 fF/EL1 ; fbL2Test = x2 fF/EL2 ; 
          bL1bL1Test = x1 wL1 fF/EL1^2 ; 
          bL2bL2Test = x2 wL2 fF/EL2^2 ; 
          bL1bL2Test = x2 wL12 fF/(EL1 EL2) ; 
          bL2bL1Test = x1 wL12 fF/(EL1 EL2) ] ; 
     If [ flag==False,           
          If [ fF < cpMin || fF > cpMax || 
               bL1fTest < cpMin || bL1fTest > cpMax || 
               bL2fTest < cpMin || bL2fTest > cpMax || 
               fbL1Test < cpMin || fbL1Test > cpMax || 
               fbL2Test < cpMin || fbL2Test > cpMax || 
               bL1bL1Test < cpMin || bL1bL1Test > cpMax || 
               bL2bL2Test < cpMin || bL2bL2Test > cpMax || 
               bL1bL2Test < cpMin || bL1bL2Test > cpMax || 
               bL2bL1Test < cpMin || bL2bL1Test > cpMax, 
               r1Error = r1 ; 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 r1 = ", r1Error//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." ] } ;

printWarning2::usage = "printWarning2 warns that a computational error has been detected (used when flag = True)." ; 
printWarning2 := { 
Print [ "WARNING:  A COMPUTATIONAL ERROR WAS DETECTED AT r1 = ", r1Error//N, "." ] ; 
Print [ "Refer to the first warning message for more info." ] } ;

printWarning3::usage = "printWarning3 warns that a computational error has been detected at r1 = zero (used when r1ZeroError = True)." ; 
printWarning3 := { 
Print [ "WARNING:  A COMPUTATIONAL ERROR WAS DETECTED AT r1 = 0." ] ; 
Print [ "Try increasing the value of the 'loopNumber' parameter (in the closed" ] ; 
Print [ "cell just above), or decreasing the value of 'totalL2' (for type 1" ] ; 
Print [ "isotherms) or 'freeL2' (for type 2 isotherms)." ] } ;

ne3::usage = "ne3 gives equations for the neighbor-effect parameters for an unsaturated lattice (used for type 1 or 2 isotherms when wL1, wL2 and wL12 are not the same)." ;

ne3 := {{EL1, EL2, r1, r2} = {eL1, eL2, R1, R2} /. 
            FindRoot[{
                eL1 (F + (R1/eL1) + (R2/eL2)) == 
                  F + ((wL1 R1)/eL1) + ((wL12 R2)/eL2), 
                eL2 (F + (R1/eL1) + (R2/eL2)) == 
                  F + ((wL2 R2)/eL2) + ((wL12 R1)/eL1), 
                (R1/(eL1^2)) == 
                  FreeL1 kL1 F 
                    ((((F/(F + (R1/eL1) + (R2/eL2)))))^(nL1 - 1)), 
                (R2/(eL2^2)) == 
                  FreeL2 kL2 F 
                    ((((F/(F + (R1/eL1) + (R2/eL2)))))^(nL2 - 1)), 
                F == 1 - nL1 R1 - nL2 R2}, {eL1, EL1}, {eL2, EL2}, {R1, 
                r1}, {R2, r2}, {F, f}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 18]; 
        f = 1 - nL1 r1 - nL2 r2; x1 = (r1/f); x2 = (r2/f)};

ne3::usage = "ne3 gives equations for the neighbor-effect parameters for an unsaturated lattice (used for type 1 or 2 isotherms when wL1, wL2 and wL12 are not the same)." ; 
(* Alternate Version *)

ne3 := {{EL1, EL2, r1, r2} = {eL1, eL2, R1, R2} /. 
            FindRoot[{
                eL1 == ((F + ((wL1 R1)/eL1) + ((wL12 R2)/eL2))/(F + 
                      (R1/eL1) + (R2/eL2))), 
                eL2 == ((F + ((wL2 R2)/eL2) + ((wL12 R1)/eL1))/(F + 
                      (R1/eL1) + (R2/eL2))), 
                (R1/(eL1^2)) == 
                  FreeL1 kL1 F 
                    ((((F/(F + (R1/eL1) + (R2/eL2)))))^(nL1 - 1)), 
                (R2/(eL2^2)) == 
                  FreeL2 kL2 F 
                    ((((F/(F + (R1/eL1) + (R2/eL2)))))^(nL2 - 1)), 
                F == 1 - nL1 R1 - nL2 R2}, {eL1, EL1}, {eL2, EL2}, {R1, 
                r1}, {R2, r2}, {F, f}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 18]; 
        f = 1 - nL1 r1 - nL2 r2; x1 = (r1/f); x2 = (r2/f)};

ne4::usage = "ne4 gives the value of the neighbor-effect parameters below saturation for type 1 or 2 isotherms when wL1 = wL2 = wL12." ;

ne4 := {{EL1, EL2, r1, r2} = {eL2, eL2, R1, R2} /. 
            FindRoot[{
                eL2 F + R1 + R2 == ((eL2 F + wL1 R1 + wL1 R2)/eL2), 
                (R1/(eL1^2)) == 
                  FreeL1 kL1 F 
                    (((((F eL2)/(F eL2 + R1 + R2))))^(nL1 - 1)), 
                (R2/(eL2^2)) == 
                  FreeL2 kL2 F 
                    (((((F eL2)/(F eL2 + R1 + R2))))^(nL2 - 1)), 
                F == 1 - nL1 R1 - nL2 R2}, {eL2, EL2}, {R1, r1}, {R2, r2}, {
                F, f}, AccuracyGoal -> 7, MaxIterations -> 200, 
              WorkingPrecision -> 18]; f = 1 - nL1 r1 - nL2 r2; 
        x1 = (r1/f); x2 = (r2/f)};

ne5::usage = "ne5 gives the value of the neighbor-effect parameters when wL1 = wL2 = wL12 = 1." ;

ne5 := {{EL1, EL2, r1, r2} = {1, 1, R1, R2} /. 
            FindRoot[{
                R1 == FreeL1 kL1 F ((((F/(F + R1 + R2))))^(nL1 - 1)), 
                R2 == FreeL2 kL2 F ((((F/(F + R1 + R2))))^(nL2 - 1)), 
                F == 1 - nL1 R1 - nL2 R2}, {R1, r1}, {R2, r2}, {F, f}, 
              AccuracyGoal -> 7, MaxIterations -> 200, 
              WorkingPrecision -> 18]; f = 1 - nL1 r1 - nL2 r2; 
        x1 = (r1/f); x2 = (r2/f)};

neL2a::usage = "neL2a gives the values of EL2 and r2 at a known concentration of free ligand 2 (used only when r1 = 0 with type 1 isotherms)." ;

neL2a := {{EL2, X2, r2} = {eL2, z, R2} /. 
            FindRoot[{eL2 (eL2 + z) == eL2 + z wL2, 
                (z/(eL2^2)) == FreeL2 kL2 (((1 + (z/eL2)))^(1 - nL2)), 
                z == (R2/(1 - nL2 R2))}, {eL2, EL2}, {z, X2}, {R2, r2}, 
              AccuracyGoal -> 7, MaxIterations -> 200, 
              WorkingPrecision -> 18]};

neL2b::usage = "neL2b gives the values of EL2 and r2 at a known concentration of free ligand 2 (used only when r1 = 0 with type 2 isotherms)." ;

neL2b := {{EL2, X2} = {eL2, z} /. 
            FindRoot[{eL2 (eL2 + z) == eL2 + z wL2, 
                (z/(eL2^2)) == FreeL2 kL2 (((1 + (z/eL2)))^(1 - nL2))}, {
                eL2, EL2}, {z, X2}, AccuracyGoal -> 7, 
              MaxIterations -> 200, WorkingPrecision -> 18]; 
        r2 = (X2/(1 + X2 nL2))};

startPrint::usage = "startPrint starts the printout of the Scatchard plot." ; 
startPrint := { 
     Print [ "concen-" ] ; 
     Print [ "tration of" ] ; 
     If [ type==1, Print [ "total L1      r1          r1/L1f" ] ] ; 
     If [ type==2, Print [ "free  L1      r1          r1/L1f" ] ] ; 
     If [ r1ZeroError==False, Print [ "0              0          ", r1L1finitial//N ] ] } ;

Expressions for conditional probabilities when r1 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." ; 
bL1f::usage = "bL1f is the probability that a randomly chosen bL1 bound ligand end is bordered by a free lattice residue." ; 
bL2f::usage = "bL2f is the probability that a randomly chosen bL2 bound ligand end is bordered by a free lattice residue." ; 
fbL1::usage = "fbL1 is the probability that a randomly chosen free lattice residue is bordered on a given side by a bL1 end of a bound ligand." ; 
fbL2::usage = "fbL2 is the probability that a randomly chosen free lattice residue is bordered on a given side by a bL2 end of a bound ligand." ; 
bL1bL1::usage = "bL1bL1 is the probability that a randomly chosen bL1 end of a bound ligand is bordered by another bL1 bound ligand end." ; 
bL2bL2::usage = "bL2bL2 is the probability that a randomly chosen bL2 end of a bound ligand is bordered by another bL2 bound ligand end." ; 
bL1bL2::usage = "bL1bL2 is the probability that a randomly chosen bL1 end of a bound ligand is bordered by a bL2 bound ligand end." ; 
bL2bL1::usage = "bL2bL1 is the probability that a randomly chosen bL2 end of a bound ligand is bordered by a bL1 bound ligand end." ;

ff := 1 / (1 + x1/EL1 + x2/EL2) ; 
bL1f [ h_] := h[[8]]/ h[[2]] /;  h[[1]] > 0 ; (* = ff/EL1 *) 
bL1f [ h_] := BL1F /;  h[[1]]==0 ; 
bL2f [ h_] := h[[8]]/ h[[3]] /;  h[[1]] > 0 ; (* = ff/EL2 *) 
bL2f [ h_] := BL2F /;  h[[1]]==0 ; 
fbL1 [ h_] := h[[6]] h[[8]]/ h[[2]] /;  h[[1]] > 0 ; (* = x1 ff/EL1 *) 
fbL1 [ h_] := FBL1 /;  h[[1]]==0 ; 
fbL2 [ h_] := h[[7]] h[[8]]/ h[[3]] /;  h[[1]] > 0 ; (* = x2 ff/EL2 *) 
fbL2 [ h_] := FBL2 /;  h[[1]]==0 ; 
bL1bL1 [ h_] := h[[6]] wL1 h[[8]]/ h[[2]]^2 /;  h[[1]] > 0 ; (* = x1 wL1 ff/EL1^2 *) 
bL1bL1 [ h_] := BL1BL1 /;  h[[1]]==0 ; 
bL2bL2 [ h_] := h[[7]] wL2 h[[8]]/ h[[3]]^2 /;  h[[1]] > 0 ; (* = x2 wL2 ff/EL2^2 *) 
bL2bL2 [ h_] := BL2BL2 /;  h[[1]]==0 ; 
bL1bL2 [ h_] := h[[7]] wL12 h[[8]]/(h[[2]] h[[3]]) /;  h[[1]] > 0 ; 
     (* = x2 wL12 ff/(EL1 EL2) *) 
bL1bL2 [ h_] := BL1BL2 /;  h[[1]]==0 ; 
bL2bL1 [ h_] := h[[6]] wL12 h[[8]]/(h[[2]] h[[3]]) /;  h[[1]] > 0 ; 
     (* = x1 wL12 ff/(EL1 EL2) *) 
bL2bL1 [ h_] := BL2BL1 /;  h[[1]]==0 ;

Define plotting routines

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

makePlot2::usage = "makePlot2 plots the theoretical r2 versus ligand 2 concentration curve (in black) along with corresponding experimental data points (in gray)." ; 
makePlot2 [ yOrigin_, xLabel_, xRange_, yRange_ ] := { 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, 
          If [ r1ZeroError, printWarning2, printWarning1 ] ] ; 
     plotTable = Table [ {abscissa[[i]], ordinate[[i,4]]}, 
          {i,points+1-dropPoints} ] ; 
     If [ r1ZeroError==False,
          curvePlot = ListPlot [ plotTable, PlotJoined->True, 
               DisplayFunction->Identity, 
               If [ defaultOutput==False, 
                    PlotRange->{{0, xRange}, {0, yRange}} , 
                    PlotRange->All ] ] ; 
          pointPlot = ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
               DisplayFunction->Identity, 
               If [ defaultOutput==False, 
                    PlotRange->{{0, xRange}, {0, yRange}} , 
                    PlotRange->All ] ] ; 
          plotTable2 = 
               Table [ {concArray[[i]], experimentalR2[[i]]}, 
                    {i,points} ]//N; 
          curvePlot2 = ListPlot [ plotTable2, PlotJoined->True, 
               AxesLabel->{xLabel, "r2"}, 
               PlotLabel->"  [bound ligand 2]/[polymer]", 
               DisplayFunction->Identity, 
               PlotStyle->{GrayLevel [0.75]}, 
               If [ defaultOutput==False, 
                    PlotRange->{{0, xRange}, {0, yRange}} , 
                    PlotRange->All ] ] ; 
          pointPlot2 = ListPlot [ plotTable2, 
               PlotStyle->{PointSize[0.015], GrayLevel [0.75] }, 
               DisplayFunction->Identity, 
               If [ defaultOutput==False, 
                    PlotRange->{{0, xRange}, {0, yRange}} , 
                    PlotRange->All ] ] ; 
          Show [ curvePlot2, pointPlot2, curvePlot, pointPlot, 
               DisplayFunction->$DisplayFunction, 
               AxesOrigin->{0,yOrigin} ] ] ; 
     N [ plotTable ] }

makeResidualsPlot::usage = "makeResidualsPlot is a plotting routine for the residuals in the plot of experimental versus theoretical r2 vs. [ligand 2] data points.  Ligand 2 concentration is plotted on the abscissa." ; 
makeResidualsPlot [ dataList_, curveList_, xLabel_ ] := { 
     residualArray = dataList-curveList ; 
     If [ saveOutput==False, 
          Unprotect [ Out ] ; Clear [ Out ] ; Protect [ Out ] ] ; 
     If [ flag, printWarning2 ] ; 
     If [ expZero==1, 
          concList = abscissa , 
          concList = Drop [ abscissa, 1 ] ] ; 
     plotTable = Table [ {concList[[j]], residualArray[[j]]}, 
          {j,points+expZero-dropPoints} ] ; 
     If [ r1ZeroError==False, 
          curvePlot := ListPlot [ plotTable, PlotJoined->True, 
               AxesLabel->{xLabel, "residuals"}, 
               PlotLabel->"exp.-theoretical", 
               DisplayFunction->Identity, 
               PlotRange->{{0, maxConcL1}, All} ] ; 
          pointPlot := ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
               DisplayFunction->Identity, 
               PlotRange->{{0, maxConcL1}, 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 r2 vs. [ligand 2] data points.  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 ] ; 
     If [ expZero==1, 
          pointList = Range [ 0, points ] , 
          pointList = Range [ points ] ] ; 
     plotTable = Table [ {pointList[[j]], residualArray[[j]]}, 
          {j,points+expZero-dropPoints} ] ; 
     If [ r1ZeroError==False, 
          curvePlot := ListPlot [ plotTable, PlotJoined->True, 
               AxesLabel->{"pt. #", "residuals"}, 
               PlotLabel->"exp.-theoretical", 
               DisplayFunction->Identity, 
               PlotRange->All ] ; 
          pointPlot := ListPlot [ plotTable, PlotStyle->PointSize[0.015], 
               DisplayFunction->Identity, 
               PlotRange->All ] ; 
          Show [ curvePlot, pointPlot, 
               DisplayFunction->$DisplayFunction] ] ; 
     N [ plotTable ] }

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

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