(* ::Package:: *)

(* Defining plotting functions *)

<<"ErrorBarPlots`";
<<"PlotLegends`"

DataPlot[xlist_, ylist_, range_, title_, axeslabel_] :=
(
g1=ListPlot[Transpose[{xlist,ylist}],PlotMarkers->{Graphics[Rectangle[{0., 0.}], ImageSize -> {11., Automatic}],0.03}];
Show[g1,PlotRange-> range,Frame->True,FrameLabel->axeslabel,
BaseStyle->{FontFamily-> "Helvetica",FontSize->11,FontWeight->"Bold"},
Axes->False,AspectRatio->0.7`,DisplayFunction->Identity,PlotLabel->title]
)

DataFunctionPlot[xlist_, ylist_,function_, range_, title_, axeslabel_] :=
(
Block[{g1,g2},
g1=Plot[function,{x,range[[1,1]],range[[1,2]]},PlotStyle->{Red,Thick}];
g2=ListPlot[Transpose[{xlist,ylist}],PlotMarkers->{Graphics[Rectangle[{0., 0.}], ImageSize -> {11., Automatic}],0.02}];
Show[g1,g2,PlotRange-> range,Frame->True,FrameLabel->axeslabel,
BaseStyle->{FontFamily-> "Helvetica",FontSize->11,FontWeight->"Bold"},
Axes->False,AspectRatio->0.7`,DisplayFunction->Identity,PlotLabel->title ]
]
)

DataFunctionPlotWithErrors[xlist_, ylist_, errlist_, function_, range_, title_, axeslabel_] :=
(
Block[{g1,g2},
g1=Plot[function,{x,range[[1,1]],range[[1,2]]},PlotStyle->{Red,Thick}];
g2=ErrorListPlot[Table[{{xlist[[j]],ylist[[j]]},ErrorBar[errlist[[j]]]},{j,1,Length[xlist]}],PlotMarkers->{Graphics[Rectangle[{0., 0.}], ImageSize -> {11., Automatic}],0.02}];
Show[g1,g2,PlotRange-> range,Frame->True,FrameLabel->axeslabel,
BaseStyle->{FontFamily-> "Helvetica",FontSize->11,FontWeight->"Bold"},
Axes->False,AspectRatio->0.7`,DisplayFunction->Identity,PlotLabel->title ]
]
)

DataListLinePlot[xlist_, ylist_, modelxylist_, range_, title_, axeslabel_] :=
(
Block[{g1,g2},
g1=ListLinePlot[modelxylist,InterpolationOrder->1,PlotStyle->{Red,Thick}];
g2=ListPlot[Transpose[{xlist,ylist}],PlotMarkers->{Graphics[Rectangle[{0., 0.}], ImageSize -> {11., Automatic}],0.02}];
Show[g1,g2,PlotRange-> range,Frame->True,FrameLabel->axeslabel,
BaseStyle->{FontFamily-> "Helvetica",FontSize->11,FontWeight->"Bold"},
Axes->False,AspectRatio->0.7`,DisplayFunction->Identity,PlotLabel->title ]
]
)

(* Defining expressions for physical quantites *)

FractionalOccupancies[Stot_,Dtot_,K1_,alpha_]:= Module[{T0roots,T1roots,T2roots,T0,T1,T2,fractions},
T0roots = NSolve[4 K1^2 Dtot^2 (1 - alpha) T0^3 + 4 K1^2 Dtot (-Dtot + 2 Dtot alpha + Stot - Stot alpha) T0^2 + (4 K1^2 Dtot Stot alpha - K1^2 Stot^2 alpha - 4 K1^2 Dtot^2 alpha - 2 K1 Stot - 1) T0 + 1 == 0];
T1roots = NSolve[K1^2 Dtot^2 alpha (1 - alpha) T1^3 + 2 K1 Dtot (1 - alpha - K1 Dtot alpha) T1^2 +(2 K1^2 Dtot Stot alpha - K1^2 Stot^2 alpha - 2 K1 Dtot - 2 K1 Stot-1) T1 + 2 K1 Stot == 0];
T2roots = NSolve[4 K1^2 Dtot^2 alpha (1 - alpha) T2^3 + 4 K1 Dtot (-2 + 2 alpha - K1 Stot alpha + K1 Dtot alpha^2 + K1 Stot alpha^2) T2^2 + alpha (4 K1^2 Dtot Stot - 4 K1^2 Dtot Stot alpha - 4 K1^2 Dtot^2 - K1^2 Stot^2 alpha - 2 K1 Stot - 4 K1 Dtot - 1) T2 + K1^2 Stot^2 alpha^2 == 0];

fractions={0.3,0.4,0.3};
Do[test={T0,T1,T2}/.Flatten[{T0roots[[i]],T1roots[[j]],T2roots[[k]]}];
If[Im[test]=={0,0,0},
If[Abs[(test[[1]]+test[[2]]+test[[3]]-1.0)]<= 0.000000001 && 1 >= test[[1]] >= 0 && 1 >= test[[2]] >= 0 && 1 >= test[[3]] >= 0,fractions=test;]
]
,{i,Length[T0roots]},{j,1,Length[T1roots]},{k,Length[T2roots]}];
Return[fractions];
]

(* Fitting algorithm - LevenbergMarquardt alternating with parameter random modification *)

DimerBindingFit[t0obs_, errort0_, t1obs_, errort1_,
        t2obs_, errort2_, ligand_, protein_, consta_] :=
      Module[{dp, dadp, dbdp, dcdp, \[Alpha], a, \[Beta], weightinga,
          weightingb, weightingc, ba, bb,
          bc, \[CapitalDelta]ia, \[CapitalDelta]ib, \[CapitalDelta]ic,
\[Chi]old, \[Chi], \[Lambda], p, \[CapitalDelta]p,
          loop, \[CapitalDelta]itesta, \[CapitalDelta]itestb,
\[CapitalDelta]itestc},

p={(Max[ligand]-Min[ligand])/2,1.0};
If[consta!=0, p[[2]]=consta];

dp = 0.00000005;
dadp = Table[0.0, {2}, {Length[ligand]}];
dbdp = Table[0.0, {2}, {Length[ligand]}];
dcdp = Table[0.0, {2}, {Length[ligand]}];
\[Alpha] = Table[0.0, {2}, {2}];
\[Beta] = Table[0.0, {2}];
If[((errort0 . errort0) >
                0.0000001) && ((errort1 . errort1) >
                0.0000001) && ((errort2 . errort2) >
                0.0000001),
weightinga = 1.0/(errort0*errort0);
weightingb = 1.0/(errort1*errort1);
weightingc = 1.0/(errort2*errort2); ,
weightinga = Table[1.0, {Length[ligand]}];
weightingb = Table[1.0, {Length[ligand]}];
weightingc = Table[1.0, {Length[ligand]}]; ];

computedvalues = Transpose[Table[
              FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]], p[[2]]], {k, 1,Length[ligand]}]];
\[CapitalDelta]ia = t0obs - computedvalues[[1]];
\[CapitalDelta]ib = t1obs - computedvalues[[2]];
\[CapitalDelta]ic = t2obs - computedvalues[[3]];

bestp = p;
best\[Chi] =  (weightinga*\[CapitalDelta]ia) . \[CapitalDelta]ia +
              (weightingb*\[CapitalDelta]ib) . \[CapitalDelta]ib +
              (weightingc*\[CapitalDelta]ic) . \[CapitalDelta]ic;

Do[
           p[[1]] = p[[1]]*RandomReal[{1.05, 0.95}];
           p[[2]] = p[[2]]*RandomReal[{1.05, 0.95}]; 
           If[consta!=0, p[[2]]=consta];

           computedvalues = Transpose[Table[
              FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]], p[[2]]], {k, 1,Length[ligand]}]];
           \[CapitalDelta]ia = t0obs - computedvalues[[1]];
           \[CapitalDelta]ib = t1obs - computedvalues[[2]];
           \[CapitalDelta]ic = t2obs - computedvalues[[3]];

            \[Chi]old = (weightinga*\[CapitalDelta]ia) .  \[CapitalDelta]ia +
                        (weightingb*\[CapitalDelta]ib) .  \[CapitalDelta]ib +
                        (weightingc*\[CapitalDelta]ic) .  \[CapitalDelta]ic;

            computedderivates = Transpose[
              (Table[FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]]+dp, p[[2]]], {k, 1,Length[ligand]}]
              - Table[FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]]-dp, p[[2]]], {k, 1,Length[ligand]}])/(2*dp)];
            dadp[[1]] = computedderivates[[1]];
            dbdp[[1]] = computedderivates[[2]];
            dcdp[[1]] = computedderivates[[3]];
            computedderivates = Transpose[
              (Table[FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]], p[[2]]+dp], {k, 1,Length[ligand]}]
              - Table[FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]], p[[2]]-dp], {k, 1,Length[ligand]}])/(2*dp)];
            dadp[[2]] = computedderivates[[1]];
            dbdp[[2]] = computedderivates[[2]];
            dcdp[[2]] = computedderivates[[3]];
  
              \[Alpha] = \[Alpha]*0.0;
              Do[Do[
                  Do[\[Alpha][[i, j]] +=
                  weightinga[[k]]*dadp[[i, k]]*dadp[[j, k]]; , {k, 1,Length[ligand]}];
                  Do[\[Alpha][[i,j]] +=
                  weightingb[[k]]*dbdp[[i, k]]*dbdp[[j, k]]; , {k, 1,Length[ligand]}];
                  Do[\[Alpha][[i,j]] += 
                  weightingc[[k]]*dcdp[[i, k]]* dcdp[[j, k]]; , {k, 1,Length[ligand]}]; 
              , {j, 1,2}]; , {i, 1,2}];

              loop =0;
              \[Lambda] = 0.1;
              \[CapitalDelta]p ={0.1,0.1,0.1};
              While[(Apply[Plus,Abs[\[CapitalDelta]p]] > 0.00001) && loop < 101, loop += 1;
              ba = dadp . (weightinga* \[CapitalDelta]ia);
              bb = dbdp . (weightingb* \[CapitalDelta]ib);
              bc = dcdp . (weightingc* \[CapitalDelta]ic);
              \[Beta] = ba + bb + bc;
              a = \[Alpha];

              Do[a[[i,i]] = (1 + \[Lambda])*\[Alpha][[i,i]], {i, 1,2}];
              \[CapitalDelta]p = Inverse[a] . \[Beta];
              If[consta!=0, \[CapitalDelta]p[[2]]=0.0];
              If[p[[1]]+\[CapitalDelta]p[[1]] < 0.0000001, \[CapitalDelta]p[[1]] = 0.0000001 - p[[1]]];
              If[p[[2]]+\[CapitalDelta]p[[2]] < 0.0000001, \[CapitalDelta]p[[2]] = 0.0000001 - p[[2]]];

              computedvalues = Transpose[Table[
              FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]]+ \[CapitalDelta]p[[1]], p[[2]]+ \[CapitalDelta]p[[2]]], {k, 1,Length[ligand]}]];
              \[CapitalDelta]itesta = t0obs - computedvalues[[1]];
              \[CapitalDelta]itestb = t1obs - computedvalues[[2]];
              \[CapitalDelta]itestc = t2obs - computedvalues[[3]];

              \[Chi] =(weightinga*\[CapitalDelta]itesta) .  \[CapitalDelta]itesta + (weightingb*\[CapitalDelta]itestb) .  \[CapitalDelta]itestb + (weightingc*\[CapitalDelta]itestc) .  \[CapitalDelta]itestc;
              If[0.999999999999 > (\[Chi]/\[Chi]old), 
                 \[Lambda] *= 0.1; \[Chi]old = \[Chi]; \[CapitalDelta]ia = \[CapitalDelta]itesta; \[CapitalDelta]ib = \[CapitalDelta]itestc; \[CapitalDelta]ic = \[CapitalDelta]itestc; p += \[CapitalDelta]p; ];
              If[1.0 < (\[Chi]/\[Chi]old), \[Lambda] *= 10.0; ]; ];
              If[\[Chi]old < best\[Chi], best\[Chi] = \[Chi]old; bestp = p; ]; 

, {20}];

 (* Compute a rough error estimate for parameters using residuals as error estimates for individual data points*)

          If[consta!=0, p[[2]]=consta];
              computedvalues = Transpose[Table[
              FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]], p[[2]]], {k, 1,Length[ligand]}]];
              \[CapitalDelta]itesta = t0obs - computedvalues[[1]];
              \[CapitalDelta]itestb = t1obs - computedvalues[[2]];
              \[CapitalDelta]itestc = t2obs - computedvalues[[3]];
              weightinga = 1.0/(\[CapitalDelta]itesta*\[CapitalDelta]itesta);
              weightingb = 1.0/(\[CapitalDelta]itestb*\[CapitalDelta]itestb);
              weightingc = 1.0/(\[CapitalDelta]itestc*\[CapitalDelta]itestc); 

            computedderivates = Transpose[
              (Table[FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]]+dp, p[[2]]], {k, 1,Length[ligand]}]
              - Table[FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]]-dp, p[[2]]], {k, 1,Length[ligand]}])/(2*dp)];
            dadp[[1]] = computedderivates[[1]];
            dbdp[[1]] = computedderivates[[2]];
            dcdp[[1]] = computedderivates[[3]];
            computedderivates = Transpose[
              (Table[FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]], p[[2]]+dp], {k, 1,Length[ligand]}]
              - Table[FractionalOccupancies[ligand[[k]],protein[[k]],p[[1]], p[[2]]-dp], {k, 1,Length[ligand]}])/(2*dp)];
            dadp[[2]] = computedderivates[[1]];
            dbdp[[2]] = computedderivates[[2]];
            dcdp[[2]] = computedderivates[[3]];
  
              \[Alpha] = \[Alpha]*0.0;
              Do[Do[
                  Do[\[Alpha][[i, j]] +=
                  weightinga[[k]]*dadp[[i, k]]*dadp[[j, k]]; , {k, 1,Length[ligand]}];
                  Do[\[Alpha][[i,j]] +=
                  weightingb[[k]]*dbdp[[i, k]]*dbdp[[j, k]]; , {k, 1,Length[ligand]}];
                  Do[\[Alpha][[i,j]] += 
                  weightingc[[k]]*dcdp[[i, k]]* dcdp[[j, k]]; , {k, 1,Length[ligand]}]; 
              , {j, 1,2}]; , {i, 1,2}];

Return[{bestp, best\[Chi],Inverse[\[Alpha]]}];];

