Basic Operating Procedure

1) Evaluate the "Initialize" notebook for the case you're interested in; then close it and open the corresponding "Output" notebook. 
2) Set the "Input Parameters" to the desired values. 
3) Evaluate the "Output" notebook. 
4) Repeat 2) and 3) as desired.

Since there is some redundancy in the various output parameters that are plotted and tabulated, you may want to unlock and make inactive or delete cells for output you're not interested in.

Introduction

These Mathematica notebooks deal with a simple model for the non-sequence-selective binding of small molecules (ligands) such as proteins to long linear polymers such as DNA.  They calculate and plot data for the binding of ligands to an infinite linear lattice for the following cases:  
     a) symmetric ligands of one type ("Sym");  
     b) symmetric ligands of two types binding to an isotropic lattice ("2Lig");  
     c) asymmetric ligands of one type binding to an isotropic lattice ("Iso");  
     d) asymmetric ligands of one type binding to an anisotropic lattice ("Ani").  
In the last case ("Ani") two binding modes are allowed, i.e., different numbers of lattice residues may be covered when the ligand binds in different orientations.  Case a is equivalent to the case analyzed by McGhee & von Hippel in their well-known paper (J. Mol. Biol. vol. 86, 469-489 [1974].  Erratum vol. 103, 679 [1976]).  Cases b-d are more general.

For each case there is an initialization notebook that contains the relevant sets of equations, and an output notebook that calculates and plots the Scatchard plot, the neighbor-effect parameters, and the conditional probabilities, all as functions of lattice saturation (r).  In addition, for all cases except the two ligand case ("2Lig"), the output notebooks calculate the 2D (at a given r) and 3D (for the whole range of lattice saturation) cluster length distributions.  The method used is based on the treatment given in Wolfe & Meehan, J. Mol. Biol. vol. 223, 1063-1087 (1992).  It is important to note that the programs for 2 types of symmetric ligand on an isotropic lattice ("2Lig") or asymmetric ligands on isotopic ("Iso") or anisotropic lattices ("Ani") will not give correct output under all circumstances.  IT IS THE RESPONSIBILITY OF THE USER TO EVALUATE THE CORRECTNESS OF HIS OR HER OUTPUT.  This can be done in a straightforward manner (see below).

Using the Programs

To use the notebooks, first open the initialization notebook for the desired case and evaluate the "initialization cells" either via the dialog box or the "Action" menu.  This needs to be done only once per Mathematica session.  Then, close the initialization notebook, open the corresponding output notebook, and specify the desired "input variables" values.  Integer values for n are expected but are generally not required.  [However, the programs for asymmetric ligands on an anisotropic lattice ("Ani") with two binding modes (n1 not equal to n2), and for type 2 or type 3 isotherms for two kinds of symmetric ligands on an isotropic lattice ("2Lig") will not work for non-integral values of n.]  The variable "points" gives the number of points at equally spaced degrees of lattice occupancy above r=0 to be used in the plot.  Because parameter values frequently change rapidly when the lattice is nearly saturated, I have added several "extra" points just before rSat (the maximum v!
 alue of r at lattice saturation), using gray text in the "initialize output arrays" cell.

The output notebook can be executed either in its entirety (with the "Evaluate Notebook" command from the "Action" menu) or in sections (by selecting the cells and hitting the "enter" key).  If the program is being executed in sections, the Scatchard plot must generally be calculated and plotted first (this can be done by selecting and executing the group of cells in the right-most bracket starting at the "input variables" cell).  The only exception to this rule is in the case of the calculation of cluster length distributions for symmetric ligands of one type ("Sym"), which can be done without any preliminary calculation of a Scatchard plot; here, it is only necessary to execute the cells grouped immediately below the input variables cell.

The Scatchard plot for an asymmetric ligand with a single binding mode (n1 = n2) binding to an isotropic ("Iso") or anisotropic ("Ani") lattice is accompanied by a "reference" symmetric ligand Scatchard plot having the same intercepts and initial slope (in gray).  Since the asymmetric ligand Scatchard plot (in black) is plotted after the reference plot, the reference plot will be invisible if the two plots are identical or nearly so.  The reference plot can be eliminated, if desired, by setting the Boolean variable "referencePlot" to "False" (in the cell just under the input variables).

After the Scatchard plot has been calculated, one or more of the neighbor-effect parameters, conditional probabilities, or cluster-length distributions can be plotted in any order.  The only restriction is that additional calculations are needed for both the 2D and 3D cluster length distributions, and these must be performed (by the cell just below the respective heading) before the data can be plotted.  The particular degree of saturation for a 2D cluster length distribution is specified by "rCluster2D", which equals r/rSat.  2D cluster length distributions can be calculated repeatedly for different r values after the Scatchard plot has been calculated once.  The number of equally spaced values of r at which the cluster distribution will be calculated for the 3D plot is specified by "points3D".  The amount of time needed for the cluster length distribution calculation increases exponentially with the length of cluster being considered (you set the maximum length for the cal!
 culations with "maxLength2D" or "maxLength3D").  However, data is printed as it is calculated, so it won't be lost if you abort.

When the computer eventually runs out of memory, it is necessary to quit Mathematica and start over.  It is also best to quit Mathematica and start over when switching from one case to another.  When using version 1 of Mathematica, don't leave the "show page breaks" option on when executing the program, as this interferes with the Print command and causes loss of output.

Most of the cells in the programs are locked, so they would have to be unlocked to make any changes.  Also, many of the cells are closed; to view their contents, open them via the "Cells" menu (Mathematica version 1) or via "Attributes" in the "Style" menu (Mathematica version 2).  The output from the programs is presented both graphically and numerically.  The numeric output may be eliminated by putting a semicolon at the end of the line in the cell just before the corresponding graph.

A meaningless equation ("scroll = 1") has been added to the end of a number of the text cells in the output notebooks.  This allows all of a notebook's results to be displayed by automatic scrolling when the entire notebook is being evaluated (the "Action - Evaluate Notebook" command).  Otherwise, the automatic scrolling usually stops at an inactive cell (at least in the Macintosh version of Mathematica).

Recognizing Incorrect Output

Problems may occasionally arise with execution of these programs, signaled by an unexpected discontinuity or plateau in one or more of the plots, a message that FindRoot has encountered a single Jacobian, or a message that "Newton's method has failed to converge to the prescribed accuracy" after the specified iteration limit.  The discontinuity or plateau may indicate that FindRoot has jumped to an incorrect solution of the equations being solved for.  This type of problem may not be evident in the Scatchard plot if it occurs at a point where the curve is close to the x-axis, but it will generally be obvious in the plots of one or more of the neighbor-effect parameters or q (particularly in the latter case).  These problems usually arise when the values of one or more of the neighbor-effect parameters (or q) are changing rapidly as a function of r, and the step size used for r in the iteration is too large (i.e., the value of "points" is too small).

If meaningless or physically impossible values are generated for one or more of the neighbor-effect parameters, q, or conditional probabilities, a flag will be set and warning messages will accompany the output.  In addition, data which has been determined to be incorrect will not be displayed.  HOWEVER, THE ABSENCE OF A WARNING MESSAGE IS NO GUARANTEE THAT THE OUTPUT IS CORRECT.  Incorrect values may still be in the physically possible range.  If the results of any calculation are in doubt, a simple and reliable way to test the correctness of the output is to change (e.g. double or quadruple) the number of points in the plot.  If the initial output was incorrect, the numbers above the r value where the mistake occurred will change in the second calculation.

Avoiding Incorrect Output

The problem of FindRoot jumping to the wrong root can generally be overcome by one of three methods:  1) by increasing the number of points in the plot; 2) by choosing less extreme values of the input variables; 3) by changing the form of the equations being solved.  The third method should be used when the first doesn't work and you don't want to use the second.  To use this approach you need to note which set of equations are being solved when the problem arises (given by the value of "ne").  Then go back into the initialization notebook and find the equations in question.  In many cases there will be a standard version of these equations, given in an initialization cell, and one or more "alternate" versions, in cells which are not initialization cells.  The version of the equations in the initialization cell is the one that normally executes most rapidly; however, one of the slower versions may be less susceptible to jumping to the wrong root in a particular situation.  S!
 elect and execute a cell containing one of the alternate versions and retry the calculation.  This approach cannot be used for errors that occur at saturation, because alternate forms of the equations used for this calculation are not available.

An example of a situation where calculating the correct output is difficult is the anisotropic lattice case with k1 = k2 = n1 = n2 = 1, w11 = 10, w22 = 0, w12 = 0.01, w21 = 0.011.  The ligands in this case initially form tightly bound head-to-head pairs, but must rearrange into a much weaker tail-to-head arrangement close to saturation.  Mistakes are often made in this calculation without setting a flag.  The error is obvious when inspecting the plot for q (and the plots for some of the b-b conditional probabilities).  The second alternative set of ne1 equations works much better for this case than the primary version.

Program Bugs and Limitations

A discontinuity that appears between the last and next-to-last points in a plot that does not respond to the above-noted measures may represent a program bug (most likely in the anisotropic lattice case ("Ani") with n1 not equal to n2).  If you believe you have encountered a bug, please report it to the programmer (Alan R. Wolfe, Dept. of Biopharmaceutical Sciences, University of California, San Francisco, CA 94143-0446; e-mail address:  alwolfe@bigfoot.com).  Remember that the parameter values will not always be correctly calculated at saturation if you use nonintegral values of n.

=======================
End of Read Me
=======================