Chemiosmotic Rate Law

Standard Form - Isolated Eepp Complex

Rate Law Definition

State variables

var={
Dre (* Reduced Donor EC Specie Concentration *),
Dtot(* Total Donor Reactant Concentration = Dre+Dox *),
Aox(* Oxidized Acceptor EC Specie Concentration *),
Atot (* Total Acceptor Reactant Concentration =Are+Aox *)
};

Constants

T=300;(* Degree Kelvin *)
R=8.31;(*(J/(K mol))*)
F=96.485;(*(C/(1000 mol))*)
n = 2 (* Number of Electrons Transfered Per Turn *);
e=0 (* Electron Crossing Within Complex *);
ΔEm=359.49 (* Standard Midpoint Electrical Potential (mV) from complex I*);
TOL=0.0001; (* Small Number *)
νdr=1;
νdo=1;
νar=1;
νao=1;
oneColumn=233.553;
twoColumns=242.03;
colorSelect=39;

Initial Condition

initial={
Dtot[0]==100,
Atot[0]==100,
Dre[0]==100-TOL,
Aox[0]==100-TOL
};

Adjustable Parameters

adjust={
Vmaxf-> 21 (* Maximum forward reaction velocity ((μM)/Min)*),
KsD -> 5 (* Dissociation Equilibrium Constant for D μM *),
KsA->17(* Dissociation Equilibrium Constant for D μM *),
KsDA->0(* Menten Constant for DA μM *),
KsAD ->0(* Menten Constant for AD μM *),
KrD->0.25(* Redox State Constant for D *),
KrA->0.018(* Redox State Constant for A *)};

Constraints

Keq=Exp[((n F ΔEm )/(R T))];
Vmaxr=(Vmaxf/Keq);
Dox[t_]:=Dtot[t]-Dre[t];
Are[t_]:=Atot[t]-Aox[t];
Dratio[t_]:=(((Dre[t])^(νdr))/((Dox[t])^(νdo)));
Aratio[t_]:=(((Aox[t])^(νao))/((Are[t])^(νar)));
 ΔΨ[t_]:=0; (* No electric potential across membrane *)

Vmaxr/.adjust (* Displays the calculated value for Vmaxr *)
Keq (* Displays the value for Keq *)

Thermodynamic Force

ΔXred[t_]:=ΔEm + ((R T)/(n F)) Log[Dratio[t]Aratio[t]]+e ΔΨ[t];
ΔX[t_]:=n F ΔXred[t];

Thermodynamic Force Function

ThermoF[t_]:=1-Exp[((-ΔX[t])/(R T))];
ThermoR[t_]:=Exp[((ΔX[t])/(R T))]-1;
(*Thermo[t_]:=If[ΔX[t]<0,(Exp[((ΔX[t])/(R T))]-1),(1-Exp[((-ΔX[t])/(R T))])];*)

Redox State Function

RedoxF[t_]:=(((Dratio[t])/(Dratio[t]+KrD)))(((Aratio[t])/(Aratio[t]+KrA)));
RedoxR[t_]:=(((1/Dratio[t])/(1/Dratio[t]+1/KrD)))(((1/Aratio[t])/(1/Aratio[t]+1/KrA)));

Saturation Function

So[t_]:=((Dtot[t]Atot[t])/(Dtot[t]Atot[t]+KsDA Dtot[t]+KsD KsDA)); (* Ordered Sequential bi bi Binding Mechanism *)
Sr[t_]:=(Dtot[t]Atot[t])/(Dtot[t]Atot[t]+KsDA Dtot[t]+KsAD Atot[t]+KsD KsDA);
(* Random Sequential bi bi Binding Mechanism *)
Sp[t_]:=(Dtot[t]Atot[t])/(Dtot[t]Atot[t]+KsD Dtot[t]+KsA Atot[t]+KsD KsA);
(* Parallel bi bi Binding Mechanism *)

Rate Expression

J[t_]:=Vmaxf Sp[t]RedoxF[t]ThermoF[t]
(*J[t_]:=If[ΔX[t]<0,Vmaxr Sp[t]Redox[t]Thermo[t],Vmaxf Sp[t]Redox[t]Thermo[t]]*)
Jr[t_]:=Vmaxr Sp[t] RedoxR[t]ThermoR[t]

Rate Law Testing

myODE={
Dre'[t]== -J[t],
Dtot'[t]==0,
Aox'[t]== -J[t],
Atot'[t]==0
};
sim= NDSolve[Join[myODE/.adjust,initial],var,{t,0,100}]

Plot[{Dox[t]/.sim,Dre[t]/.sim},{t,0,20},AxesLabel->{"Time (min)","Velocity (((μM)/Min))"}]

Plot[{Aox[t]/.sim,Are[t]/.sim},{t,0,20},AxesLabel->{"Time (min)","Velocity (((μM)/Min))"}]

Plot[{J[t]/.sim/.adjust,Jr[t]/.sim/.adjust},{t,0,20}]

Loading Experimental Data

File Input

SetDirectory[NotebookDirectory[]];

SetDirectory["ComplexI_121207"]

filenames1=FileNames["1-*.txt"]
filenames2=FileNames["2-*.txt"]
filenames3=FileNames["3-*.txt"]
filenames4=FileNames["4-*.txt"]
filenames5=FileNames["5-*.txt"]
filenames6=FileNames["6-*.txt"]
filenames7=FileNames["7-*.txt"]
filenames8=FileNames["8-*.txt"]
filenames9=FileNames["9-*.txt"]
filenames10=FileNames["10-*.txt"]
filenames11=FileNames["11-*.txt"]
filenames12=FileNames["12-*.txt"]
filenames13=FileNames["13-*.txt"]
filenames14=FileNames["14-*.txt"]
filenames15=FileNames["15-*.txt"]
filenames16=FileNames["16-*.txt"]

Donor Data Array

data100d=Take[ImportString[StringReplace[Import[filenames1[[1]]] ,","->"."],"Table"], {11,-2}];
data75d=Take[ImportString[StringReplace[Import[filenames7[[1]]] ,","->"."],"Table"], {11,-2}];
data50d=Take[ImportString[StringReplace[Import[filenames6[[1]]] ,","->"."],"Table"], {11,-2}];
data25d=Take[ImportString[StringReplace[Import[filenames3[[1]]] ,","->"."],"Table"], {11,-2}];
data10d=Take[ImportString[StringReplace[Import[filenames4[[1]]] ,","->"."],"Table"], {11,-2}];
data5d=Take[ImportString[StringReplace[Import[filenames5[[1]]] ,","->"."],"Table"], {11,-2}];

Acceptor Data Array

data150a=Take[ImportString[StringReplace[Import[filenames16[[1]]] ,","->"."],"Table"], {11,-2}];
data100a=Take[ImportString[StringReplace[Import[filenames11[[1]]] ,","->"."],"Table"], {11,-2}];
data75a=Take[ImportString[StringReplace[Import[filenames15[[1]]] ,","->"."],"Table"], {11,-2}];
data50a=Take[ImportString[StringReplace[Import[filenames12[[1]]] ,","->"."],"Table"], {11,-2}];
data25a=Take[ImportString[StringReplace[Import[filenames13[[2]]] ,","->"."],"Table"], {11,-2}];
data12a=Take[ImportString[StringReplace[Import[filenames14[[2]]] ,","->"."],"Table"], {11,-2}];

Data Preprocessing

Acceptor Curves

Definitions

(rawDataFunctionA[ToExpression[#]]=Interpolation[Take[SetAccuracy[ToExpression[StringJoin["data",#,"a"]],5],{1,-1,1}],InterpolationOrder->10,Method->"Spline"])&/@{"150","100","75","50","25","12"};

(endPointA[ToExpression[#]]=IntegerPart[rawDataFunctionA[ToExpression[#]][[1,1]][[2]]-10])&/@{"150","100","75","50","25","12"};

(endTimeA[ToExpression[#]]=IntegerPart[rawDataFunctionA[ToExpression[#]][[1,1]][[2]]/60.0])&/@{"150","100","75","50","25","12"};

Finding the Slope of the Background Reaction and Correcting the Raw Data Function

(slopeA[ToExpression[#]]=(1/((IntegerPart[endPointA[ToExpression[#]] 0.7])))(Abs[rawDataFunctionA[ToExpression[#]][endPointA[ToExpression[#]]]-rawDataFunctionA[ToExpression[#]][endPointA[ToExpression[#]]-(IntegerPart[endPointA[ToExpression[#]] 0.7])]]))&/@{"75","50"};

(slopeA[ToExpression[#]]=(1/((IntegerPart[endPointA[ToExpression[#]] 0.7])))(Abs[rawDataFunctionA[ToExpression[#]][endPointA[ToExpression[#]]]-rawDataFunctionA[ToExpression[#]][endPointA[ToExpression[#]]-(IntegerPart[endPointA[ToExpression[#]] 0.7])]]))&/@{"25","12"};

(slopeA[ToExpression[#]]=(1/((IntegerPart[endPointA[ToExpression[#]] 0.2])))(Abs[rawDataFunctionA[ToExpression[#]][endPointA[ToExpression[#]]]-rawDataFunctionA[ToExpression[#]][endPointA[ToExpression[#]]-(IntegerPart[endPointA[ToExpression[#]] 0.2])]]))&/@{"150","100"};

(*slopeA[100]=slopeA[50]*)
slopeA[150]=slopeA[75];
slopeA[50]=slopeA[75];
slopeA[25]=slopeA[75];
slopeA[12]=slopeA[75];

(correctedRawDataFunctionA[ToExpression[#]][p_]=rawDataFunctionA[ToExpression[#]][p]+(slopeA[ToExpression[#]] p))&/@{"150","100","75","50","25","12"};

Finding the offset

(offSetA[ToExpression[#]]=SetAccuracy[Mean[Commonest[Take[Table[correctedRawDataFunctionA[ToExpression[#]][p],{p,0,endPointA[ToExpression[#]]}],{IntegerPart[(endPointA[ToExpression[#]] -0.1)],-1,1}],3]],4])&/@{"150","100","75","50","25","12"};

Finding the offset (Beta)

(offSetA[ToExpression[#]]=SetAccuracy[Min[Table[correctedRawDataFunctionA[ToExpression[#]][p],{p,0,endPointA[ToExpression[#]]}]],4])&/@{"150","100","75","50","25","12"};

Finding the Vertical Shift from 100uM (Depreciated)

(startShiftA[ToExpression[#]]=correctedRawDataFunctionA[100][0]-correctedRawDataFunctionA[ToExpression[#]][0])&/@{"150","100","75","50","25","12"};

(shiftA[ToExpression[#]]=correctedRawDataFunctionA[100][endPointA[100]]-correctedRawDataFunctionA[ToExpression[#]][endPointA[ToExpression[#]]])&/@{"150","100","75","50","25","12"};

Shifting the Data; Conversion from Absorbtion to [Concentration];  Correcting the time scale

endTimeA[12]=5;
endTimeA[25]=6;
endTimeA[50]=10;
endTimeA[75]=12;

(concDataFunctionA[ToExpression[#]][t_]:=(correctedRawDataFunctionA[ToExpression[#]][t*60]-offSetA[ToExpression[#]])*1000/6.22/;t<=endTimeA[ToExpression[#]])&/@{"150","100","75","50","25","12"};

(*(tempConcDataFunctionA[ToExpression[#]][t_]:=(correctedRawDataFunctionA[ToExpression[#]][t*60]-offSetA[ToExpression[#]])*1000/6.22/;t<=endTimeA[ToExpression[#]])&/@{"150"};*)

(*concDataFunctionA[150][t_]:=tempConcDataFunctionA[150][t]*(((concDataFunctionA[100][0])/(tempConcDataFunctionA[150][0])))*)

Normalized Curves

(normalConcDataFunctionA[ToExpression[#]][t_]:=concDataFunctionA[ToExpression[#]][t]/concDataFunctionA[ToExpression[#]][0])&/@{"150","100","75","50","25","12"};

Normalized Curves for CoQ Simulation

(normalCoQConcDataFunctionA[ToExpression[#]][t_]:=1-normalConcDataFunctionA[ToExpression[#]][t])&/@{"150","100","75","50","25","12"};

Donor Curves

Definitions

(rawDataFunction[ToExpression[#]]=Interpolation[Take[SetAccuracy[ToExpression[StringJoin["data",#,"d"]],5],{1,-1,1}],InterpolationOrder->10,Method->"Spline"])&/@{"100","75","50","25","10","5"};

(endPoint[ToExpression[#]]=IntegerPart[rawDataFunction[ToExpression[#]][[1,1]][[2]]-10])&/@{"100","75","50","25","10","5"}

(endTime[ToExpression[#]]=IntegerPart[rawDataFunction[ToExpression[#]][[1,1]][[2]]/60.0])&/@{"100","75","50"}
endTime[25]=5;
endTime[10]=3;
endTime[5]=3;

Finding the Slope of the Background Reaction and Correcting the Raw Data Function

(slope[ToExpression[#]]=(1/((IntegerPart[endPoint[ToExpression[#]] 0.3])))(Abs[rawDataFunction[ToExpression[#]][endPoint[ToExpression[#]]]-rawDataFunction[ToExpression[#]][endPoint[ToExpression[#]]-(IntegerPart[endPoint[ToExpression[#]] 0.3])]]))&/@{"100","75","50","25","10","5"}

slope[100]=slopeA[100];
(*slope[100]=slopeA[50]
slope[75]=slopeA[50]
slope[50]=slopeA[50]
slope[25]=slopeA[50]
slope[10]=slopeA[50]
slope[5]=slopeA[50]*)

(correctedRawDataFunction[ToExpression[#]][p_]=rawDataFunction[ToExpression[#]][p]+(slope[ToExpression[#]] p))&/@{"100","75","50","25","10","5"};

Finding the offset

(offSet[ToExpression[#]]=SetAccuracy[Mean[Commonest[Take[Table[correctedRawDataFunction[ToExpression[#]][p],{p,0,endPoint[ToExpression[#]]}],{IntegerPart[(endPoint[ToExpression[#]] -0.1)],-1,1}],3]],4])&/@{"100","75","50","25","10","5"}

(*(offSet[ToExpression[#]]=SetAccuracy[Min[Table[correctedRawDataFunction[ToExpression[#]][p],{p,0,endPoint[ToExpression[#]]}]],4])&/@{"100","75","50","25","10","5"}*)

Shifting the curves; Conversion from Absorbtion to [Concentration];  Correcting the time scale

(concDataFunction[ToExpression[#]][t_]:=(correctedRawDataFunction[ToExpression[#]][t*60]-offSet[ToExpression[#]])*1000/6.22/;t<endTime[ToExpression[#]])&/@{"100","75","50","25","10","5"};
(*concDataFunction[100][t_]:=(correctedRawDataFunction[100][t*60]-end)*1000/6.22;*)

Normalized Potential Curves

(normalConcDataFunction[ToExpression[#]][t_]:=concDataFunction[ToExpression[#]][t]/concDataFunction[ToExpression[#]][0])&/@{"100","75","50","25","10","5"};

Vmaxf and Saturation Function Parameters Fitting

Instantenous Velocities

V100d=(concDataFunction[100][0]-concDataFunction[100][1]);
V75d=(concDataFunction[75][0]-concDataFunction[75][1]);
V50d=(concDataFunction[50][0]-concDataFunction[50][1]);
V25d=(concDataFunction[25][0]-concDataFunction[25][1]);
V10d=(concDataFunction[10][0]-concDataFunction[10][1]);
V5d=(concDataFunction[5][0]-concDataFunction[5][1]);

V100a=(concDataFunctionA[100][0]-concDataFunctionA[100][1]);
V75a=(concDataFunctionA[75][0]-concDataFunctionA[75][1]);
V50a=(concDataFunctionA[50][0]-concDataFunctionA[50][1]);
V25a=(concDataFunctionA[25][0]-concDataFunctionA[25][1]);
V12a=(concDataFunctionA[12][0]-concDataFunctionA[12][1]);

Linear Fit of Starting Concentration

linearFit=Flatten[{FindFit[Table[{t,concDataFunction[100][t]},{t,0.4,1.4,0.1}], -V100 t+ b100, {V100,b100},t],FindFit[Table[{t,concDataFunction[75][t]},{t,0,1,0.1}], -V75 t+ b75, {V75,b75},t],FindFit[Table[{t,concDataFunction[50][t]},{t,0.2,endTime[50]*0.075,0.1}], -V50 t+ b50, {V50,b50},t],FindFit[Table[{t,concDataFunction[25][t]},{t,0,endTime[25]*0.075,0.1}], -V25 t+ b25, {V25,b25},t],FindFit[Table[{t,concDataFunction[10][t]},{t,0,endTime[10]*0.075,0.1}], -V10 t+ b10, {V10,b10},t],FindFit[Table[{t,concDataFunction[5][t]},{t,0,endTime[5]*0.075,0.1}], -V5t+ b5, {V5,b5},t]
}];
Show[{Plot[{(-V100 t+b100)/.linearFit,(-V75 t+ b75)/.linearFit,(-V50 t+ b50)/.linearFit,(-V25 t+ b25)/.linearFit,(-V10 t+ b10)/.linearFit},{t,0,15},PlotStyle->{{Green,Dashed}},PlotRange->{{0,10},{0,100}},AxesLabel->{"Time (Min)", "NADH (μM)"},BaseStyle -> {Thick,FontFamily->"Helvetica", FontWeight -> "Bold", FontSize -> 12},ImageSize -> Large]},{Plot[{concDataFunction[100][t],concDataFunction[75][t],concDataFunction[50][t],concDataFunction[25][t],concDataFunction[10][t]},{t,0,15},PlotStyle->ColorData[colorSelect,"ColorList"],AxesLabel->{"Time (Min)", "NADH (μM)"},BaseStyle -> {Thick,FontFamily->"Helvetica", FontWeight -> "Bold", FontSize -> 12},ImageSize -> Large]},PlotRange->{{0,10},{0,100}}]

linearFitA=Flatten[{FindFit[Table[{t,concDataFunctionA[150][t]},{t,0,endTimeA[150]*0.05,0.1}], -V150A t+ b150A, {V150A,b150A},t],FindFit[Table[{t,concDataFunctionA[100][t]},{t,0,endTimeA[100]*0.05,0.1}], -V100A t+ b100A, {V100A,b100A},t],FindFit[Table[{t,concDataFunctionA[75][t]},{t,0,endTimeA[75]*0.05,0.1}], -V75A t+ b75A, {V75A,b75A},t],FindFit[Table[{t,concDataFunctionA[50][t]},{t,0,endTimeA[50]*0.05,0.1}], -V50A t+ b50A, {V50A,b50A},t],FindFit[Table[{t,concDataFunctionA[25][t]},{t,0,endTimeA[25]*0.05,0.1}], -V25A t+ b25A, {V25A,b25A},t],FindFit[Table[{t,concDataFunctionA[12][t]},{t,0,endTimeA[12]*0.05,0.1}], -V12A t+ b12A, {V12A,b12A},t]
}];
Show[Plot[{(-V100A t+b100A)/.linearFitA,(-V75A t+ b75A)/.linearFitA,(-V50A t+ b50A)/.linearFitA,(-V25A t+ b25A)/.linearFitA},{t,0,15},PlotStyle->{{Green,Dashed}},BaseStyle -> {Thick,FontFamily->"Helvetica", FontWeight -> "Bold", FontSize -> 12},ImageSize -> Large],Plot[{concDataFunctionA[100][t],concDataFunctionA[75][t],concDataFunctionA[50][t],concDataFunctionA[25][t],concDataFunctionA[12][t]},{t,0,15},PlotRange->{-0.5,100},LabelStyle->{Medium,(FontFamily->"Helvetica")},AxesLabel->{"Time (Min)", "ΔNADH (μM)"},PlotStyle->{ColorData[colorSelect,"ColorList"][[7]],ColorData[colorSelect,"ColorList"][[2]],ColorData[colorSelect,"ColorList"][[3]],ColorData[colorSelect,"ColorList"][[4]],ColorData[colorSelect,"ColorList"][[5]],ColorData[colorSelect,"ColorList"][[6]]},BaseStyle -> {Thick,FontFamily->"Helvetica", FontWeight -> "Bold", FontSize -> 12},ImageSize -> Large],PlotRange->{{0,15},{0,100}}]

Donor Saturation Parameter (Subscript[K, S])

dataD={{b100A,V100A},{b100,V100},{b75,V75}, {b50,V50}, {b25,V25}, {b10,V10}}/.linearFit/.linearFitA;
param=FindFit[dataD,((VmaxD S)/(KsD+S)),{VmaxD,KsD},{S}];
nlmD=NonlinearModelFit[dataD,((VmaxD S)/(KsD+S)),{VmaxD,KsD},{S}];
nlmD[{"BestFit","RSquared","ParameterTable"}]
StdD=nlmD["ParameterTable"][[1,1]][[2,3]];

Donor Saturation Parameter (Subscript[K, S])

dataA={{b100A,V100A},{b100,V100}, {b75A,V75A},{b50A,V50A}, {b25A,V25A},{b12A,V12A}}/.linearFit/.linearFitA
paramA=FindFit[dataA,((VmaxA S)/(KsA+S)),{VmaxA,KsA},{S}];
nlmA=NonlinearModelFit[dataA,((VmaxA S)/(KsA+S)),{VmaxA,KsA},{S}];
nlmA[{"BestFit","RSquared","ParameterTable"}]

StdA=nlmA["ParameterTable"][[1,1]][[2,3]];

Combined Vmaxf

Vmaxf=((VmaxD+VmaxA)/2)/.param/.paramA
PooringStandardError=Sqrt[(((VmaxD-1)(StdD^2)+(VmaxA-1)(StdA^2))/((VmaxD-1)+(VmaxA-1)))]/.param/.paramA

Michaelis - Menten Plot

fittedPlot=Plot[{((VmaxD S)/(KsD+S))/.param,Vmaxf}, {S, 0,100},PlotStyle->{Darker[Blue],{Darker[Cyan],Thick,Dashed}},BaseStyle -> {FontFamily->"Arial", FontSize -> 8, PrivateFontOptions->{"OperatorSubstitution"->False}},
PlotRange->{{0,101},{0,22}},Frame->{{True,False},{True,False}},FrameLabel->{"NADH (μM)","Velocity (μM/Min)"}];
dataPlot=ListPlot[dataD,
PlotRange->{{0,100},{0,22}},
PlotStyle->Red,BaseStyle -> {FontFamily->"Arial", FontSize -> 8, PrivateFontOptions->{"OperatorSubstitution"->False}},
PlotMarkers->{Graphics[{Red,Text[StyleForm["o",FontFamily->"Arial",FontSize->12,FontWeight->"Bold"]]}]}
];


michaelismenten=Show[{fittedPlot,dataPlot,Graphics[Text[StyleForm["(Subscript[V, max])f",FontFamily->"Arial",FontSize->10,FontWeight->"Bold"],{12,20},{0,0},{1,0}]],Graphics[Text[StyleForm["(Subscript[K, S, D])",FontFamily->"Arial",FontSize->10,FontWeight->"Bold"],{(KsD/.param)-1.5,2},{-1,0},{0,1}]],Graphics[{Dashed,Thick,Purple,Line[{{KsD/.param,0},{KsD/.param,(VmaxD/2.)/.param}}]}]}]


fittedPlotA=Plot[{((VmaxA S)/(KsA+S))/.paramA,Vmaxf}, {S, 0,151},PlotStyle->{Darker[Blue],{Darker[Cyan],Thick,Dashed}},BaseStyle -> {FontFamily->"Arial", FontSize -> 8, PrivateFontOptions->{"OperatorSubstitution"->True,"WindowsUseTrueTypeNames"->True}},
PlotRange->{{0,101},{0,22}},Frame->{{True,False},{True,False}},FrameLabel->{"CoQ (μM)","Velocity (μM/Min)"}];
dataPlotA=ListPlot[dataA,
PlotRange->{{0,100},{0,22}},
PlotStyle->Orange,BaseStyle -> {FontFamily->"Arial", FontSize -> 8, PrivateFontOptions->{"OperatorSubstitution"->False,"WindowsUseTrueTypeNames"->True}},
PlotMarkers->{Graphics[{Orange,Text[StyleForm["o",FontFamily->"Arial",FontSize->12,FontWeight->"Bold"]]}]}];
michaelismentenA=Show[{fittedPlotA,dataPlotA,Graphics[Text[StyleForm["(Subscript[V, max])f",FontFamily->"Arial",FontSize->10,FontWeight->"Bold"],{18,20},{0,0},{1,0}]],Graphics[Text[StyleForm["(Subscript[K, S, A])",FontFamily->"Arial",FontSize->10,FontWeight->"Bold"],{(KsA/.paramA)-2.5,2},{-1,0},{0,1}]],Graphics[{Dashed,Thick,Purple,Line[{{KsA/.paramA,0},{KsA/.paramA,(VmaxA/2.)/.paramA}}]}]}]

Redox State Function Parameter Fitting

Supporting Functions

model[V_?NumberQ,ksd_?NumberQ,ksa_?NumberQ,krd_?NumberQ,kra_?NumberQ,gamma_?NumberQ,Dtotal_?NumberQ,Dneg_?NumberQ,Atotal_?NumberQ,Aplus_?NumberQ,runTime_?NumberQ]:=(model[V,ksd,ksa,krd,kra,gamma,Dtot,Dneg,Atot,Aplus]=First[Dre/.NDSolve[Join[{
Dre'[t]== -J[t],
Dtot'[t]==0,
Aox'[t]== -J[t],
Atot'[t]==0
} /.{
Vmaxf->V,KsD->ksd,KsA->ksa,KrD->krd,KrA->kra,γ->gamma
},{
Dtot[0]==Dtotal,
Atot[0]==Atotal,
Dre[0]==Dneg,
Aox[0]==Aplus
}], {Dtot,Atot,Dre,Aox}, {t,runTime}]]);

modelA[V_?NumberQ,ksd_?NumberQ,ksa_?NumberQ,krd_?NumberQ,kra_?NumberQ,gamma_?NumberQ,Dtotal_?NumberQ,Dneg_?NumberQ,Atotal_?NumberQ,Aplus_?NumberQ,runTime_?NumberQ]:=(modelA[V,ksd,ksa,krd,kra,gamma,Dtot,Dneg,Atot,Aplus]=First[Aox/.NDSolve[Join[{
Dre'[t]== -J[t],
Dtot'[t]==0,
Aox'[t]== -J[t],
Atot'[t]==0
} /.{
Vmaxf->V,KsD->ksd,KsA->ksa,KrD->krd,KrA->kra,γ->gamma
},{
Dtot[0]==Dtotal,
Atot[0]==Atotal,
Dre[0]==Dneg,
Aox[0]==Aplus
}], {Dtot,Atot,Dre,Aox}, {t,runTime}]]);

test[krd_,kra_,Atotal_]:=model[V100/.linearFit,0,0,krd,kra,1,(b100+TOL)/.linearFit,b100/.linearFit,Atotal,Atotal-TOL,endTime[100]]
test75[krd_,kra_,Atotal_]:=model[V75/.linearFit,0,0,krd,kra,1,(b75+TOL)/.linearFit,b75/.linearFit,Atotal,Atotal-TOL,endTime[75]]
test50[krd_,kra_,Atotal_]:=model[V50/.linearFit,0,0,krd,kra,1,(b50+TOL)/.linearFit,b50/.linearFit,Atotal,Atotal-TOL,endTime[50]]
test25[krd_,kra_,Atotal_]:=model[V25/.linearFit,0,0,krd,kra,1,(b25+TOL)/.linearFit,b25/.linearFit,Atotal,Atotal-TOL,endTime[25]]

concDataPlots=Plot[{concDataFunction[100][t],concDataFunction[75][t],concDataFunction[50][t],concDataFunction[25][t]},{t,0,12}];

Manipulate[Show[{Plot[Evaluate[test[krd,kra,Aox][t]],{t,0,20},PlotStyle->{Green}],Plot[concDataFunction[100][t],{t,0,20}]},PlotRange->{0,100}],{krd,0,1,0.01},{kra,0,1,0.01},{Aox,100-TOL,0,1}]

Fitting Result

multiSimulationDataD[time_]:={model[V100/.linearFit,0,0,krd,kra,1,(b100+TOL)/.linearFit,b100/.linearFit,100,100-TOL,time*endTime[100]][time*endTime[100]],0.5model[V75/.linearFit,0,0,krd,kra,1,(b75+TOL)/.linearFit,b75/.linearFit,100,100-TOL,(time*endTime[75])][time*endTime[75]],0.5model[V50/.linearFit,0,0,krd,kra,1,(b50+TOL)/.linearFit,b50/.linearFit,100,100-TOL,time*endTime[50]][time*endTime[50]],0.5model[V25/.linearFit,0,0,krd,kra,1,(b25+TOL)/.linearFit,b25/.linearFit,100,100-TOL,time*endTime[50]][time*endTime[50]]}
multiExperimentalDataD[time_]={concDataFunction[100][time*endTime[100]],0.5concDataFunction[75][time*endTime[75]],0.5concDataFunction[50][time*endTime[50]],concDataFunction[50][time*endTime[50]]};
multiErr[t_]:=Total[(multiSimulationDataD[t]-multiExperimentalDataD[t])^2]
multiNewdat=Table[{t,0},{t,0.4,0.9,0.1}];
(*Fitting the 0.4 to 0.9 range where the redox state parameters have the most effect *)

(* Parameter Fitting Routine. Uncomment to run *)
(*nlmMultiFitResult=NonlinearModelFit[multiNewdat,{multiErr[t], {0<=krd<1,0<=kra<1}},{krd,kra},t];
param=nlmMultiFitResult["BestFitParameters"];*)

param={krd->0.304,kra->0.05}; (*saved results to avoid recomputation everytime*)

(* Show Plot Results of the Fit *)
Manipulate[Show[{Plot[Evaluate[test[krd,kra,Aox][t]],{t,0,20},PlotStyle->{Green},PlotRange->{0,100}],Plot[Evaluate[test75[krd,kra,Aox][t]],{t,0,10},PlotStyle->{Green},PlotRange->{0,100}],Plot[Evaluate[test50[krd,kra,Aox][t]],{t,0,10},PlotStyle->{Green},PlotRange->{0,100}],Plot[Evaluate[test25[krd,kra,Aox][t]],{t,0,5},PlotStyle->{Green},PlotRange->{0,100}],concDataPlots},LabelStyle->{8,(FontFamily->"Arial")},Frame->{{True,False},{True,False}},FrameLabel->{"Time (Min)", "[NADH] (μM)"}],{{krd,krd/.param,"KRD Variation"},0,1,0.01},{{kra,kra/.param,"KRA Variation"},0,1,0.01},{{Aox,100,"A total Concentration"},100,0,1}]