The Binding of Two Kinds of Symmetric Ligands to an Infinite Isotropic Linear Lattice - Exp. vs. Theoretical Data by Alan R. Wolfe

This Mathematica  notebook calculates and plots the Scatchard plot, neighbor-effect parameters, conditional probabilities, and other types of data for the binding of two kinds of symmetric ligands to an isotropic lattice.  It can also compare experimental data with the theoretical binding curve and compute a correlation coefficient for the two.  It can only be used after the appropriate initialization notebook has been evaluated.  To perform the calculations, simply set the input variables to the desired values and execute the "Kernel-> Evaluation-> Evaluate Notebook" command.  This may be repeated with different values for the input variables.  

The notebook may also be executed in sections.  To operate in this manner, the input parameters cell (as well as the experimental data cell, if desired) must be executed first, followed by the cell immediately below the "Calculate…" heading.  Any of the plots can then be displayed in any order.  It's easiest to start by simply selecting the right-most cell grouping bracket and executing everything from the input parameters to the r2 versus ligand 1 concentration plot.  

The calculations are performed over a range of ligand 1 concentrations (not for an array of r values of ligand 1), and do not reach lattice saturation.  The experimental data needs to be given as an array of concentrations of ligand 1 (either total concentrations or free concentrations, depending on whether type is set to 1 or 2, respectively), and an array of r values for ligand 2 (bound ligand concentration divided by the polymer concentration in terms of monomer units) for some constant concentration of either total or free ligand 2 (again depending on the value of type).  The initial value of experimentalR2 must be positive, and experimentalData must be True, in order to perform the experimental vs. theoretical comparison.  If experimentalData is False, only theoretical data from zero up to [ligand 1] = maxConcL1 will be generated (over n = points equal intervals).  If experimentalData is True, and the initial value of experimentalR2 is negative, again only theoretical data will be generated, but this time using the list of experimental ligand 1 concentrations provided (concArray).

INPUT VARIABLES:

kL1 = 4.9; kL2 = 1000; nL1 = 4; nL2 = 8; 
wL1 = 1; wL2 = 1; wL12 = 1; type = 1; 
defaultOutput = True ; xRange = 10 ; yRange = 0.105 ; 
(* type 1 isotherms only *)  totalL2 = 0.24; polymer = 2.4; 
(* type 2 isotherms only *)  freeL2 = 2; 
(* relevant to whether experimental data is used *) 
points = 12; maxConcL1 = 12; experimentalData = True;

"type" indicates which type of isotherm is desired.  
type = 1 -> a constant total concentration of ligand 2 and varying known concentrations of total ligand 1; 
type = 2 -> a constant concentration of free ligand 2 and varying known concentrations of free ligand 1.  
          
totalL2 is the total concentration of ligand 2; polymer is the total lattice residue concentration.  Both concentrations are expressed in the same units used for kL2. 

freeL2 is the concentration of free ligand 2 expressed in the same units used for kL2.  

xRange and yRange are only used for the first plot, the plot of r2 vs. concentration of ligand 1, when experimental and theoretical data are being compared, and defaultOutput = False.  If defaultOutput = True, all the data points will be displayed.  

Set experimentalData to True if a set of experimental ligand 1 concentrations and r2 data points are to be compared with theoretical data for the same ligand 1 concentrations.  Otherwise, a set of equally spaced concentrations (the number being determined by the value of "points") up to maxConcL1 will be used, for theoretical data only.  If experimentalData = True and the initial experimental r2 value is negative, the experimental data will be ignored, but the theoretical data will still be calculated for the given set of experimental ligand 1 concentrations (= concArray).

Experimental Data (* disregarded if experimentalData is False *)

(* The number of points in the two arrays must be equal. *)

(* sample data - experimental r values for ligand 2.  This data will be disregarded if the initial value is negative *) 

experimentalR2 = {0.0981, 0.0973, 0.0977, 0.0979, 0.0984, 0.0983, 0.0985, 0.0982, 0.0985, 0.0992, 0.0993, 0.0992, 0.099, 0.0985, 0.0978, 0.0962, 0.0948, 0.0926, 0.0904, 0.0881, 0.0851, 0.0826, 0.079, 0.0755, 0.0725, 0.0688, 0.0644, 0.0601, 0.0556, 0.051, 0.0455, 0.0385, 0.0307, 0.0238, 0.0187, 0.0144, 0.0108, 0.00835, 0.00733, 0.00511, 0.00578, 0.00435, 0.00255, 0.00266, 0.00273} ;

(* Calculate the NE parameters and Scatchard plot *)

(* Note - if the first (r1 = 0) data point on the plots is calculated incorrectly , it may help to increase the value of the variable loopNumber in the cell below this one.  Using the other version of the 2 types of symmetric ligand/isotropic lattice program (the version that does not compare experimental and theoretical data), this data point may be checked by calculating a type 3 isotherm at the calculated r2 of the initial point of the type 1 or 2 isotherm; if correct, the initial neighbor-effect values will be the same.  *)

If [ experimentalR2[[1]]<0, useData = False, useData = True ] ; 

(* Initialize program parameters and calculate binding of ligand 2 at r1 = 0 *)  
flag = False ; 
freeL1 = 0 ; r1 = 0 ; EL1 = 1 ; 
If [ wL1 == 0, r1Sat = 1/(nL1 + 1), r1Sat = 1/nL1 ] ; 
If [ wL2 == 0, r2Sat = 1/(nL2 + 1), r2Sat = 1/nL2 ] ; 
setLimits ; 
If [ type==1, FreeL1 := totalL1-R1 polymer ; 
      FreeL2 := TotalL2-R2 polymer ] ; 
If [ type==2, FreeL1 := freeL1 ] ; 
loopNumber = 40 ; 
If [ type==1, loopArray = Range [1, loopNumber] ; 
     X2 = 0.01 ; r2 = 0.01 ; EL2 = 1 ;  plotLabel = "[T L1]" ; 
     Map [ Function [ zz, 
          TotalL2 = totalL2 zz/loopNumber ; neL2a ], loopArray ]  ] ; 
If [ type==2, loopArray = Range [1, loopNumber] ; 
     X2 = 0.01 ; EL2 = 1 ; plotLabel = "[F L1]" ; 
     Map [ Function [ zz, 
          FreeL2 = freeL2 zz/loopNumber ; neL2b ], loopArray ]  ] ; 
r2initial = r2 ; 
f = 1 - nL1 r1 - nL2 r2 ; 
x1 = r1/f ; x2 = r2/f ; 
If [ !wL12==1, EL1 = ( EL2 + x2 wL12 )/( EL2 + x2 ) ] ; 
checkLimits ; 
If [ flag, printWarning3 ; r1ZeroError = True; dropPoints = 1, 
     r1ZeroError = False ; dropPoints = 0 ] ; 
 
(* Initialize loop parameters and output arrays *) 
neChoose ; 
r1L1finitial = r1L1f ; EL1initial = EL1 ; EL2initial = EL2 ; 
r2L2finitial = r2L2f ; 
x2initial = x2 ; 
ffinitial = ff ; 
latticeCoverageInitial = latticeCoverage ; 
If [ experimentalData==False, 
     xAxisInterval = maxConcL1/(points) ; 
     abscissa = xAxisInterval Range [1, points] ] ; 
If [ experimentalData==True, 
     If [ concArray[[1]]==0, 
          abscissa = Drop [ concArray, 1 ] ; expZero = 1 , 
          abscissa = concArray ; expZero = 0 ] ; 
     points = Length [abscissa] ] ; 
     
(* Evaluate neighbor-effect parameters for r1 < r1Sat *)
startPrint ; 
ordinate = Map [ Function [ zz, 
     If [ type==1, totalL1 = zz ] ; 
     If [ type==2, freeL1 = zz ] ; 
     f = 1 - nL1 r1 - nL2 r2 ; 
     x1 = r1/f ; 
     x2 = r2/f ; 
     If [ needCheck, 
          If [ flag==True, {EL1,EL2} = {1,1} ] ; 
          If [ flag==False, ne ; checkLimits ] ; 
          If [ flag==True, dropPoints = dropPoints + 1 ] , 
          ne ] ; 
     If [ flag==False, Print [ zz//N, "     ", r1//N, "       ", r1L1f//N ] ] ; 
     {r1L1f, EL1, EL2, r2, r2L2f,x1,x2,ff,r1,latticeCoverage} ], 
     abscissa ] ; 
     
(* Add output parameters for r1 = 0 *)
maxConcL1 = abscissa[[points]] ; 
PrependTo [ abscissa, 0 ] ; 
PrependTo [ ordinate, 
     {r1L1finitial,EL1initial,EL2initial,r2initial, 
          r2L2finitial,0,x2initial,ffinitial,0,latticeCoverageInitial} ] ;

(* Indicate which equations were used *)  scroll = 1

Print [ "ne = ", neChoice] ; 
Print [ "r1Sat = ", r1Sat//N, "; r2Sat = ", r2Sat//N ] ; 
Print [ "maximum concentration of L1 = ", maxConcL1//N ] ; 
Print [ "r1 at maximum concentration of L1 = ", r1//N ] ; 
Print [ "r2 in absence of L1 = ", r2initial//N ] ; 
Print [ "r2 at maximum concentration of L1 = ", r2//N ] ;

(* Plot r2 versus concentration of ligand 1 *)  scroll = 1

If [ experimentalData==True && useData==True, 
     makePlot2 [ r2initial, plotLabel, xRange, yRange ] ,
     makePlot [ 4, r2initial, plotLabel, "r2", 
          "  [bound ligand 2]/[polymer]" ] ]

If [ experimentalData==True && useData==True, 
     theoreticalR2 = Take [ Table [ ordinate[[i,4]], {i,points+1} ], 
          -points - expZero ] ; 
     Print [ "correlation coefficient (R) = ", 
          corrCoeff [ experimentalR2, theoreticalR2 ] ] ; 
     makeResidualsPlot [ experimentalR2, theoreticalR2, plotLabel ] ]

If [ experimentalData==True && useData==True, 
     makeResidualsPlot2 [ experimentalR2, theoreticalR2 ] ] ;

(* Plot r1/[free ligand 1] vs. conc. of ligand 1 *)  scroll = 1

makePlot [ 1, r1L1finitial, plotLabel, "r1/L1f", "        for ligand 1" ] ;

(* Plot r1 versus concentration of ligand 1 *)  scroll = 1

makePlot [ 9, 0, plotLabel, "r1", "  [bound ligand1]/[polymer]" ] ;

(* Plot the neighbor-effect parameters *)  scroll = 1

makePlot [ 2, EL1initial, plotLabel, "eL1","NE for free L1 sites" ]

makePlot [ 3, EL2initial, plotLabel, "eL2", "NE for free L2 sites" ]

(* Plot r2/[free ligand 2] versus conc. of ligand 1 *)  scroll = 1

makePlot [ 5, r2L2finitial, plotLabel, "r2/L2f","              for ligand 2" ]

(* Plot lattice coverage vs. total conc. of ligand 1 *)  scroll = 1

makePlot [ 10, latticeCoverageInitial, plotLabel, "l.c.", 
"  fraction of lattice covered" ]

(* Plot the conditional probabilities *)  scroll = 1

(* ff *)

makePlot [ 8, 0, plotLabel, "ff", "" ]

(* bL1f *)

makeCPplot [ bL1f, plotLabel, "     bL1f" ]

(* bL2f *)

makeCPplot [ bL2f, plotLabel, "     bL2f" ]

(* fbL1 *)

makeCPplot [ fbL1, plotLabel, "     fbL1" ]

(* fbL2 *)

makeCPplot [ fbL2, plotLabel, "     fbL2" ]

(* bL1bL1 *)

makeCPplot [ bL1bL1, plotLabel, "          bL1bL1" ]

(* bL2bL2 *)

makeCPplot [ bL2bL2, plotLabel, "          bL2bL2" ]

(* bL1bL2 *)

makeCPplot [ bL1bL2, plotLabel, "          bL1bL2" ]

(* bL2bL1 *)

makeCPplot [ bL2bL1, plotLabel, "          bL2bL1" ]