(* The software performs the DSR analysis of ligand-binding curves described in the following publication (submitted): "Decomposing Complex Cooperative Ligand Binding into Simple Components: Connections between Microscopic and Macroscopic Models", by Alexey Onufriev and G. Matthias Ullmann. *) (* Requirements: Mathematica 4.0 or higher *) (* Development Team: Gustavo Moura, Niteesh Bharara, and Siddharth Saraiya *) (* Supervision: Alexey Onufriev *) (* Department of Computer Science, Virginia Tech *) (* contact: onufriev@cs.vt.edu *) (* -----------------------SET UP SECTION----------------------------------- *) (* SET THE NUMBER OF BINDING SITES. *) n=4; (* READ-IN EXPERIMENTAL DATA *) (* It is assumed that the input data has the format [mu, X(mu)] *) (* mu -> ligand chemical potential (ln[activity], in units of kT), X(mu) -> aver age receptor occupancy, 0 < X < N, where N is the user-specified number of bindi ng sites. See below if input data has different format *) data = ReadList["carpHb.phosphate.Y.vs.mu.pH=6.61.OUE.fit", {Number, Number}] (*--------------------------END SET UP -------------------------------------*) Print[] Print["......working......" ] Print[] (* load the Nonlinear fit package from the Statistics Library *) < BestFitParameters], {i, 1, n/2+1}]; (*Caclulate Chi^2 for each C. *) Errors = Table[ Apply[Plus, ( Results[Exp[datax]][[i]] - datay )^2], {i,1,n/2+1}]; (*Make a final list of Chi^2, C values and parameters, and sort them in order of Chi^2 - smallest to largest*) FinalList = Sort[ Transpose[Join[{Errors}, {Cvals}, {ParamResults} ]]] (*Print the information for the C that gives best fit*) Print[] Print[] Print["---------------------------RESULTS------------------------------------"] Print[] Print[] (*Set precision *) precision = 10^-4 (* This piece of code compares Chi^2 for each C, and takes the minimum which corresponds to the correct C-value. There may be two (or more) C-values that are identical, within the precicion specified above. This likely means that some of the 2-ns order quasi-groups can be reduced to pairs of quasi-sites with positve binding constants. If that is the case, pick the smallest C. This is not very safe, though, given noisey data, and one really should check if the corresponding second-order binding polynomial has two real roots, i.e. alpha^2 - 4*gamma >= 0 *) SmallestError = FinalList[[1]][[1]]; OptimalC = FinalList[[1]][[2]]; If [ n > 1 && (FinalList[[2]][[1]] - SmallestError) < precision, If [ FinalList[[2]][[2]] < OptimalC, OptimalC = FinalList[[2]][[2]] ]; If [ n > 3 && (FinalList[[3]][[1]] - SmallestError) < precision, If [ FinalList[[3]][[2]] < OptimalC, OptimalC = FinalList[[3]][[2]] ] ] ] Print["Number of ligand-binding sites, N= ", n] Print[] Print["Best C value: ", FinalList[[1]][[2]]] Print[] Print["Chi^2 for the best C: ", FinalList[[1]][[1]]] If[ FinalList[[1]][[2]] == OptimalC-1, Print["Chi^2 for the next best C: ", FinalList[[1]][[1]]] ] If[ FinalList[[2]][[2]] == OptimalC-1, Print["Chi^2 for the next best C: ", FinalList[[2]][[1]]] ] If[ FinalList[[3]][[2]] == OptimalC-1, Print["Chi^2 for the next best C: ", FinalList[[3]][[1]]] ] Print[] Print[] Print["DSR Fit Parameters: ", FinalList[[1]][[3]]] Print[] (*compare data with the fit *) bestC = FinalList[[1]][[2]] plotfit = Plot[Results[Exp[x]][[bestC+1]], {x,-10,10}, PlotStyle-> {Hue[1], Thickness[0.005]}, DisplayFunction -> Identity]; plotdata= ListPlot[Transpose [{datax,datay}], PlotStyle->{Thickness[1]},DisplayFunction->Identity]; (*Shows the two plots from above on the same axis.*) Show[plotfit,plotdata, PlotLabel-> FontForm["red-> best DSR fit; black-> data", {"Times", 16}], Frame -> True, GridLines -> Automatic, FrameLabel -> { FontForm["chemical potential [kT]", {"Times", 16}], FontForm["receptor occupancy", {"Times", 16}] }, DefaultFont -> {"Helvetica", 14 }, DisplayFunction-> $DisplayFunction] (*Create a parameter list, with our initial fit parameters*) Plist = Table[ paramlist[[OptimalC + 1]][[i]] /. ParamResults[[OptimalC + 1]][[i]], {i, 1, 4}]; (******************************************************************************************) (*In the following we search for quasi-groups of order > 2 *) (* Need to do this if we find a negative alpha or gamma among the DSR fit parameters*) (* N=4 case. Check to see if any of the alphas are negative. If one is negative, we combine terms in the DSR decomposition until an irreducible polynomial in the denominator is obtained, with all positive coefficients corresponding to quasi-group binding sites. *) If[n == 4, If [ OptimalC == 1 , K1 = Plist[[1]]; K2 = Plist[[2]]; a = Plist[[3]]; g = Plist[[4]]; (*Check for negative parameters. If one is found, check to see what transformatoin applies*) If[a < 0, If[K1 >= Abs[a], If[g >= Abs[ a*K1], Trans13K1 = 1, Trans13K1 = 0], Trans13K1 = 0 ]; If [K2 >= Abs [a] && Trans13K1 == 0, If[g >= Abs [a*K2], Trans13K2 = 1, Trans13K2 = 0;], Trans13K2 = 0; ]; If [Trans13K1 == 0 && Trans13K2 == 0, Trans4 = 1, Trans4 = 0 ]; (* do our transformations now ... *) If[Trans13K1 == 1, paramlist = {{K2}, {K1 + a, (a*K1 + g) / (K1 + a), (K1*g) / (a*K1 + g)}}; grouplist = {2, {1, 3}}]; If[Trans13K2 == 1, paramlist = {{K1}, {K2 + a, (a*K2 + g) / (K2 + a), (K2*g) / (a*K2 + g)}}; grouplist = {2, {1, 3}}]; If[Trans4 == 1, paramlist = {{(a + K1 + K2), (g + a*K1 + a*K2 + K1*K2)/(a + K1 + K2), (g*K1 + g*K2 + a*K1*K2)/(g + a*K1 + a*K2 + K1*K2), (g*K1*K2)/(g*K1 + g*K2 + a*K1*K2)}}; grouplist = {1, {4}}], (*This will set the new parameter list, and a group list which will be used for outputting at the end*) paramlist = {{K1}, {K2}, {a, g/a}}; grouplist = {3, {1, 1, 2}}] ]; (*Next case, when the optimal value is C=2*) If [OptimalC == 2, a1 = Plist[[1]]; a2 = Plist[[2]]; g1 = Plist[[3]]; g2 = Plist[[4]]; If [ a1 < 0 || a2 < 0, Trans4 = 1, Trans4 = 0 ]; (* do our transformations now ... *) If [Trans4 == 1, paramlist = {{ (a1 + a2), (a1*a2 + g1 + g2)/(a1 + a2), (a2*g1 + a1*g2)/(a1*a2 + g1 + g2), (g1*g2)/(a2*g1 + a1*g2) }}; grouplist = {1, {4}}, paramlist = {{a1, g1/a1}, {a2, g2/a2}}; grouplist = {2, {2, 2}} ] ]; (*If the optimal C is 0, then we do have an irreducible representation *) If [OptimalC == 0, paramlist = { {Plist [[1]] }, {Plist[[2]]}, {Plist[[3]]}, {Plist[[4]]} }; grouplist = {4, {1, 1, 1, 1}} ] ]; (*This is the N=3 case*) If [n == 3, (*This is when optimal c=1 case*) If [OptimalC == 1, k = Plist[[1]]; a = Plist[[2]]; g = Plist[[3]]; (*Check to see if alpha is less than zero, if it is then do the transformation*) If [ a < 0, paramlist = { {(a + k), (g + a*k)/(a + k), (g*k)/(g + a*k)}}; grouplist = { 1, {3}}, paramlist = { {k}, {a, g/a}}; grouplist = { 2, {1, 2}} ] ]; (*This is C=0 case, no transformation is needed*) If [OptimalC == 0, paramlist = {{Plist[[1]]}, {Plist[[2]]}, {Plist[[3]]}}; grouplist = {3, {1, 1, 1}} ] ]; (*This is the N=2 case*) If [n == 2, (*This is C=1 case*) If [OptimalC == 1, paramlist = {{Plist[[1]], Plist[[2]]/Plist[[1]]}}; grouplist = {1, {2}}, paramlist = {{Plist[[1]]}, {Plist[[2]]}}; grouplist = {2, {1, 1}} ] ]; (*This is the N=1 case, just for consistency. *) If [n == 1, paramlist = {{Plist[[1]]}}; grouplist = {1, {1}} ]; (*****************************************************************************) (*Output quasi-group decomposition and constants *) Print[] Print[] Print["Quasi-group decomposition (irreducible representation): ", grouplist[[2]] ] Print[]; Print["Quasi-group binding constants, one set per group:"] (*This will be used to tell if we have aleady output a variable*) Subscript[Ksub, 1] = 1; Subscript[Ksub, 2] = 1; Subscript[Ksub, 3] = 1; Subscript[Ksub, 4] = 1; (* quasi-group binding constants. All are positive. *) For[i = 1, i <= grouplist[[1]], Print["m = ", grouplist[[2]][[i]] ]; For[k = 1, k <= grouplist[[2]][[i]], subs = ToString[ grouplist[[2]][[i]] ] <> ToString[k]; Ksubprint = Subscript[Ksub, grouplist[[2]][[i]] ]; If[subs == "11", subs = "'"]; Print[Subscript[Superscript[K, subs], Subscript[Ksub, grouplist[[2]][[i]]]], " = ", paramlist[[i]][[k]] ]; Clear[subs]; k++;]; Print[]; Print[]; Subscript[Ksub, grouplist[[2]][[i]]] = Subscript[Ksub, grouplist[[2]][[i]]] + 1; i++;]