RedoxMech - Oxidoreductase (EC 1.) xCellerator Language Extension

Examples and Guide

Package Loading

Needs["xlr8r`"](*Required for deterministic simulations*)
Needs["xSSA`"](*Required for stochastic simulations*)

<<RedoxMech`

Parameters and Constants

Molecular Species/State Variables

Donor={Dre, Dox};
Acceptor={Aox,Are};
Enz={Epp,Enp,Epn,Enn};

System Constants

R=8.31;(*(J/(K mol))*)
F=96.485;(*(C/(1000 mol))*)
nA=6.022 (10^23);(* Avocadro's number 1/mol *)
TOL=(10^(-9)); (* Small Number For Computational Purpose *)
T=300;(* Degree Kelvin *)
Volume=(10^(-18));(*1 attoL - volume of a mitochondria *)
ΔΨ[t_]:=0; (* No electric potential across membrane *)

System Initial Condition

Set1={
complexNumber->ToMolecules[30],
DreNumber->ToMolecules[2400],
DoxNumber->0,
AoxNumber->ToMolecules[2400],
AreNumber->0
};

Complex Specific Constants

Measured and Referenced Biochemical Kinetic Values

(Mitochondrial complex I - NADH dehydrogenase)

ComplexI={
kCat->0.3,(*tri-molecular catalysis*)
kCatb->0.3,(*bi-molecular catalysis*)
kCatu->0.3,(*uni-molecular catalysis*)
boundaryTransfer->0.3 (10^3),(* reaction order of the boundary electron transfer reactions *)
BindingOrder->0.3 (10^8)(* The order of the binding steps (M^(-1))*),
Ksd->12 (10^(-6))(*M*),(* Binding constant for donor *)
Ksad->12 (10^(-6))(*M*), (* Binding constant for donor after binding of acceptor *)
Ksa->16.5 (10^(-6))(*M*), (* Binding constant for acceptor *)
Ksda->16.5 (10^(-6))(*M*), (* Binding constant for acceptor after binding of donor *)
Krd->0.01,(* Tri-molecular boundary electron transfer with the donor EC *)
Krdb->0.01,(* Bi-molecular boundary electron transfer with the donor EC *)
Kra->0.1,(* Tri-molecular boundary electron transfer with the acceptor EC *)
Krab->0.1,(* Bi-molecular boundary electron transfer with the acceptor EC *)
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*)
νdr-> 1,
νdo-> 1,
νar-> 1,
νao-> 1
};

(Mitochondrial complex III)

ComplexIII={
kCat->1.0,(*tri-molecular catalysis*)
kCatb->1.0,(*bi-molecular catalysis*)
kCatu->1.0,(*uni-molecular catalysis*)
boundaryTransfer->1.0 (10^3),(* reaction order of the boundary electron transfer reactions *)
BindingOrder->1.0 (10^9)(* The order of the binding steps (M^(-1))*),
Ksd->150 (10^(-6))(*M*),(* Binding constant for donor *)
Ksad->150 (10^(-6))(*M*), (* Binding constant for donor after binding of acceptor *)
Ksa->8 (10^(-6))(*M*), (* Binding constant for acceptor *)
Ksda->8 (10^(-6))(*M*), (* Binding constant for acceptor after binding of donor *)
Krd->0.55,(* Tri-molecular boundary electron transfer with the donor EC *)
Krdb->0.55,(* Bi-molecular boundary electron transfer with the donor EC *)
Kra->2.25,(* Tri-molecular boundary electron transfer with the acceptor EC *)
Krab->2.25,(* Bi-molecular boundary electron transfer with the acceptor EC *)
n -> 2, (* Number of Electrons Transfered Per Turn *)
e-> 1,(* Electron Crossing Within Complex *)
ΔEm-> 169,(* Standard Midpoint Electrical Potential (mV) from complex III*)
νdr-> 1,
νdo-> 1,
νar-> 1,
νao-> 1
};

(Mitochondrial complex IV)

ComplexIV={
kCat->0.6,(*tri-molecular catalysis*)
kCatb->0.6,(*bi-molecular catalysis*)
kCatu->0.6,(*uni-molecular catalysis*)
boundaryTransfer->0.6 (10^3),(* reaction order of the boundary electron transfer reactions *)
BindingOrder->0.6 (10^9)(* The order of the binding steps (M^(-1))*),
Ksd->5.15 (10^(-6))(*M*),(* Binding constant for donor *)
Ksad->5.15 (10^(-6))(*M*), (* Binding constant for donor after binding of acceptor *)
Ksa->1.0 (10^(-6))(*M*), (* Binding constant for acceptor *)
Ksda->1.0 (10^(-6))(*M*), (* Binding constant for acceptor after binding of donor *)
Krd->3.56,(* Tri-molecular boundary electron transfer with the donor EC *)
Krdb->3.56,(* Bi-molecular boundary electron transfer with the donor EC *)
Kra->0.1,(* Tri-molecular boundary electron transfer with the acceptor EC *)
Krab->0.1,(* Bi-molecular boundary electron transfer with the acceptor EC *)
n -> 2, (* Number of Electrons Transfered Per Turn *)
e-> 0,(* Electron Crossing Within Complex *)
ΔEm-> 319,(* Standard Midpoint Electrical Potential (mV) from complex I*)
νdr-> 1,
νdo-> 1,
νar-> 1,
νao-> 1
};

Calculated Redox Potential and Equilibrium Values

Emrd=-Log[Krd]((R T)/(n F))(*mV*);
Emra=-Log[Kra]((R T)/(n F))(*mV*);
Keq=Exp[((n F)/(R T))(ΔEm-Emrd-Emra)];
Vmaxr=(Vmaxf/Keq);

Environmental Settings for Steady-State Macro Model

initial={
Dre[0]==ToConcentration[(DreNumber)],
Dox[0]==TOL,
Aox[0]==ToConcentration[AoxNumber],
Are[0]==TOL
};
adjust={
Vmaxf->kCat ToConcentration[complexNumber](* Maximum forward reaction velocity ((μM)/Min)*),
KsD -> Ksd (10^6)(* Dissociation Equilibrium Constant for D μM *),
KsA->Ksa (10^6)(* Dissociation Equilibrium Constant for D μM *),
KsDA->Ksda(* MM Constant for DA μM *),
KsAD ->Ksad(* MM Constant for AD μM *),
KrD->Krd(* Redox State Constant for D *),
KrA->Kra(* Redox State Constant for A *),
Dtot-> ToConcentration[DoxNumber+DreNumber],(* Total Donor Reactant Concentration = Dre+Dox *)
Atot->ToConcentration[AoxNumber+AreNumber](* Total Acceptor Reactant Concentration =Are+Aox *)
};

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

myODE={
Dre'[t]== -J[t],
Dox'[t]==J[t],
Aox'[t]== -J[t],
Are'[t]== J[t]
};

myODEt={
Dre'[t]== -Jt[t],
Dox'[t]==Jt[t],
Aox'[t]== -Jt[t],
Are'[t]== Jt[t]
};

Environmental Settings/Initial Conditions for Elementary Reaction Micro Model

microInitialConditions={Dre->DreNumber,Dox->DoxNumber,Aox->AoxNumber, Are->AreNumber, Epp->complexNumber};

General Reaction Mechanism Framework for EC 1.

Micro Model - General Set of Elementary Reactions within a Prototypical EC 1.

xCellerator Arrow Input Notation

microReactions=EC1[(Overscript[{Dre,Dox}⇄{Are,Aox}, E]),ComplexI,ReactionType->"Complete",Output->"xCellerator"];
microReactions//TableForm
(* Both Options ReactionType and Output are of default value *)

Reaction Matrix Output

To aide the visualization of the internal elementary reaction network, a n x n matrix of n unique enzyme states is constructed from the reaction list in which the sparse rate constants connects the transition from the i-th enzyme state to the j-th enzyme state (Option: Output -> “Matrix”). The rate constants can be shown in terms of the unitary kinetic rate labels, steady-state rate relations, or numerical values through the option “ShowConstants -> Unitary,” “ShowConstants -> SteadyState,” or “ShowConstants -> Values”.

completeReaction=EC1[(Overscript[{Dre,Dox}⇄{Are,Aox}, E]),ComplexI,ReactionType->"Complete",Output->"Matrix"];

ReactionMatrix[completeReaction,ShowConstants->"Unitary"]

ReactionMatrix[completeReaction,ShowConstants->"SteadyState"]

ReactionMatrix[completeReaction,ShowConstants->"Values"]

Macro Model - General Steady State Rate Equation (Standard Form of Chemiosmotic Rate Law)

Thermodynamic Force

Dratio[t_]:=(((Dre[t])^(νdr))/((Dox[t])^(νdo)));
Aratio[t_]:=(((Aox[t])^(νao))/((Are[t])^(νar)));

ΔXred[t_]:=ΔEm-Emrd-Emra + ((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;

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 Atot)/(Dtot Atot+KsDA Dtot+KsD KsDA)); (* Ordered Sequential bi bi Binding Mechanism *)
Sr[t_]:=((Dtot Atot)/(Dtot Atot+KsDA Dtot+KsAD Atot+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);*)
Sp[t_]:=((Dtot Atot)/(Dtot Atot+KsD Dtot+KsA Atot+KsD KsA));
(* Parallel bi bi Binding Mechanism *)
Spt[t_]:=(Dre[t]Aox[t])/(Dre[t]Aox[t]+KsD Dre[t]+KsA Aox[t]+KsD KsA);
Sptr[t_]:=(Dox[t]Are[t])/(Dox[t]Are[t]+KsD Dox[t]+KsA Are[t]+KsD KsA);

Rate Expression

J[t_]:=Vmaxf Sp[t]RedoxF[t]ThermoF[t]
Jt[t_]:=Vmaxf Spt[t]RedoxF[t]ThermoF[t]
Jr[t_]:=Vmaxr Sp[t] RedoxR[t]ThermoR[t]
Jtr[t_]:=Vmaxr Sptr[t] RedoxR[t]ThermoR[t]

Specific Sub-Mechanisms of EC 1.

In a general-to-specific manner, specific mechanisms of EC 1. can be obtained from the list of general internal elementary reactions. The following examples show the expanded form of common ternary enzyme kinetic mechanisms in the context of EC 1.  catalyzed reaction.

PingPong Mechanism

PingPong mechanism partitions the overall ternary enzyme reaction into: (1) an ordered binding and releasing of the first substrate and product, which leaves a modified form of the enzyme; (2) a second ordered binding and releasing of the second substrate and product, which restores the original enzyme form. Because of the introduction of the two internal reaction center in the EC 1., there are two possible variations of this mechanism: (1) the modified enzyme form undergoes internal electron transfer before complexing with the second substrate; or (2) the enzyme complexes with the second substrate before undergoing internal electron transfer.

xCellerator Arrow Input Notation

PingPongReactions=EC1[(Overscript[{Dre,Dox}⇄{Are,Aox}, E]),ComplexI,ReactionType->"PingPong"];
PingPongReactions//TableForm

Reaction Matrix Output

pingPongMatrix=EC1[(Overscript[{Dre,Dox}⇄{Are,Aox}, E]),ComplexI,ReactionType->"PingPong",Output->"Matrix"];

ReactionMatrix[pingPongMatrix,ShowConstants->"Unitary"]

ReactionMatrix[pingPongMatrix,ShowConstants->"SteadyState"]

ReactionMatrix[pingPongMatrix,ShowConstants->"Values"]

Steady-State Veocity Equation

A steady-state velocity equation specific to the sub-mechanism can be derived from the reaction matrix using the method describe by Cornish-Bowden. A modified version of Mulquiney et al’s implementation of the method is provided in RedoxMech as the function SSRateEqation.

pingPongEquationOutput=SSVelocityEq[pingPongMatrix,ShowConstants->"Unitary",Report->False]
(* Gives the velocity equation in terms of the unitary rate constants *)

pingPongValueOutput=SSVelocityEq[pingPongMatrix,ShowConstants->"Values",Report->True]
(* Replace the unitary rate constants with parameter values *)

pingPongReportButton=Button["Click to Show Steady-State Report for PingPong Mechanism",CreateDocument[SSReport/.pingPongValueOutput,WindowTitle->"PingPong Steady State Velocity Equation Report"]]

Macro Model ODE Definition from the Steady-State Velocity Equation

JC1p[t_]=flux/.pingPongValueOutput/.{Dre->Dre[t],Dox->Dox[t],Aox->Aox[t],Are->Are[t],e0->ToConcentration[complexNumber]};

myODEp={
Dre'[t]== -JC1p[t],
Dox'[t]==JC1p[t],
Aox'[t]== -JC1p[t],
Are'[t]== JC1p[t]
};

Ordered Sequential Mechanism

The ordered mechanism prescribes a specific order of substrate binding to form the central ternary catalysis complex, followed by a specific order of product unbinding from the ternary complex.

xCellerator Arrow Input Notation

OrderedSeqReactions=EC1[(Overscript[{Dre,Dox}⇄{Are,Aox}, E]),ComplexI,ReactionType->"Ordered",Output->"xCellerator"];
OrderedSeqReactions//TableForm

Reaction Matrix Output

orderedMatrix=EC1[(Overscript[{Dre,Dox}⇄{Are,Aox}, E]),ComplexI,ReactionType->"Ordered",Output->"Matrix"];

ReactionMatrix[orderedMatrix,ShowConstants->"Unitary"]

ReactionMatrix[orderedMatrix,ShowConstants->"SteadyState"]

ReactionMatrix[orderedMatrix,ShowConstants->"Values"]

Steady-State Veocity Equation

Off[Solve::ratnz];

orderedValueOutput=SSVelocityEq[orderedMatrix,ShowConstants->"Values",Report->True];

orderedReportButton=Button["Click to Show Steady-State Report for Ordered Sequential Mechanism",CreateDocument[SSReport/.orderedValueOutput,WindowTitle->"Ordered Sequential Steady State Velocity Equation Report"]]

Macro Model ODE Definition from the Steady-State Velocity Equation

JC1o[t_]=flux/.orderedValueOutput/.{Dre->Dre[t],Dox->Dox[t],Aox->Aox[t],Are->Are[t],e0->ToConcentration[complexNumber]}

myODEo={
Dre'[t]== -JC1o[t],
Dox'[t]==JC1o[t],
Aox'[t]== -JC1o[t],
Are'[t]== JC1o[t]
};

Random Sequential Mechanism

In the random sequential mechanism, the sequence order of the substrates binding and products unbinding is random, which effectively encompass all four variations of the ordered sequential mechanism simultaneously.

xCellerator Arrow Input Notation

RandomSeqReactions=EC1[(Overscript[{Dre,Dox}⇄{Are,Aox}, E]),ComplexI,ReactionType->"Random",Output->"xCellerator"];
RandomSeqReactions//TableForm

Reaction Matrix Output

randomMatrix=EC1[(Overscript[{Dre,Dox}⇄{Are,Aox}, E]),ComplexI,ReactionType->"Random",Output->"Matrix"];

ReactionMatrix[randomMatrix,ShowConstants->"Unitary"]

ReactionMatrix[randomMatrix,ShowConstants->"SteadyState"]

ReactionMatrix[randomMatrix,ShowConstants->"Values"]

Steady-State Veocity Equation

randomValueOutput=SSVelocityEq[randomMatrix,ShowConstants->"Values",Report->True];

randomReportButton=Button["Click to Show Steady-State Report for Random Sequential Mechanism",CreateDocument[SSReport/.randomValueOutput,WindowTitle->"Random Sequential Steady State Velocity Equation Report"]]

Macro Model ODE Definition from the Steady-State Velocity Equation

JC1r[t_]=flux/.randomValueOutput/.{Dre->Dre[t],Dox->Dox[t],Aox->Aox[t],Are->Are[t],e0->ToConcentration[complexNumber]};

myODEr={
Dre'[t]== -JC1r[t],
Dox'[t]==JC1r[t],
Aox'[t]== -JC1r[t],
Are'[t]== JC1r[t]
};

Deterministic Simulation Comparisons Between the Micro and Macro Forms of Each Mechanism Models

Shared Control Parameters

startTime=0;(*seconds*)
endTime=600; (*seconds*)

Micromodel Simulations

Complete Reactions

microsim=run[microReactions,{0,endTime},rates->{},initialConditions->microInitialConditions/.Set1,plot->False];

Ping-Pong

microPingPongSim=run[PingPongReactions,{0,endTime},rates->{},initialConditions->microInitialConditions/.Set1,plot->False];

Ordered Sequential

microOrderedSim=run[OrderedSeqReactions,{0,endTime},rates->{},initialConditions->microInitialConditions/.Set1,plot->False];

Random Sequential

microRandomSim=run[RandomSeqReactions,{0,endTime},rates->{},initialConditions->microInitialConditions/.Set1,plot->False];

Macromodel Simulations

Complete Reactions

macrosim= NDSolve[Join[myODE/.adjust/.ComplexI/.Set1,initial/.ComplexI/.Set1],var,{t,startTime,endTime}];

macrosimt= NDSolve[Join[myODEt/.adjust/.ComplexI/.Set1,initial/.ComplexI/.Set1],var,{t,startTime,endTime}];

Ping-Pong

macrosimp= NDSolve[Join[myODEp/.adjust/.ComplexI/.Set1,initial/.ComplexI/.Set1],var,{t,startTime,endTime}];

Ordered Sequential

macrosimo= NDSolve[Join[myODEo/.adjust/.ComplexI/.Set1,initial/.ComplexI/.Set1],var,{t,startTime,endTime}];

Random Sequential

macrosimr= NDSolve[Join[myODEr/.adjust/.ComplexI/.Set1,initial/.ComplexI/.Set1],var,{t,startTime,endTime}];

Substrate Concentration Time-Series Graphical Comparison

Free Substrate Concentration

Concentration time-series composed only of free reactant concentrations (e.g. Dre, Dox, Are, Aox) not bound to the enzyme complex

freeSubstrateCompareButton=Button["Click to Show Comparison of Micro/Macro Output-Free Reactants Concentrations Time-Series",CreateDocument[{Manipulate[Grid[{{
Plot[{
If[Complete,ToConcentration[Dox[t]]/.microsim,Null],
If[OrderedSeq,ToConcentration[Dox[t]]/.microOrderedSim,Null],
If[RandomSeq, ToConcentration[Dox[t]]/.microRandomSim,Null],
If[PingPong, ToConcentration[Dox[t]]/.microPingPongSim,Null],
If[Complete,(Dox[t])/.macrosim,Null],
If[OrderedSeq,(Dox[t])/.macrosimo,Null],
If[RandomSeq,(Dox[t])/.macrosimr,Null],
If[PingPong,(Dox[t])/.macrosimp,Null]
},{t,0,time},PlotStyle->{Directive[Green,Thick],Directive[Blue,Thick],Directive[Orange,Thick],Directive[Red,Thick],Directive[Black,Thick,Dashed],Directive[Darker[Blue],Thick,Dotted],Directive[Darker[Orange],Thick,Dotted],Directive[Darker[Red],Thick,Dotted]},Frame->{{True,False},{True,False}},FrameLabel->{Style["Time (S)",Bold], Style["Concentration (μM)",Bold]},FrameStyle->{{AbsoluteThickness[1],White},{AbsoluteThickness[1],White}},BaseStyle -> {FontFamily->"Arial", FontSize -> 12, PrivateFontOptions->{"OperatorSubstitution"->False,"WindowsUseTrueTypeNames"->True}},ImageSize->650,PlotLabel->"Dox",PlotRange->All],
Plot[{
If[Complete,ToConcentration[Dre[t]]/.microsim,Null],
If[OrderedSeq,ToConcentration[Dre[t]]/.microOrderedSim,Null],
If[RandomSeq, ToConcentration[Dre[t]]/.microRandomSim,Null],
If[PingPong, ToConcentration[Dre[t]]/.microPingPongSim,Null],
If[Complete,(Dre[t])/.macrosim,Null],
If[OrderedSeq,(Dre[t])/.macrosimo,Null],
If[RandomSeq,(Dre[t])/.macrosimr,Null],
If[PingPong,(Dre[t])/.macrosimp,Null]
},{t,0,time},PlotStyle->{Directive[Green,Thick],Directive[Blue,Thick],Directive[Orange,Thick],Directive[Red,Thick],Directive[Black,Thick,Dashed],Directive[Darker[Blue],Thick,Dotted],Directive[Darker[Orange],Thick,Dotted],Directive[Darker[Red],Thick,Dotted]},Frame->{{True,False},{True,False}},FrameLabel->{Style["Time (S)",Bold], Style["Concentration (μM)",Bold]},FrameStyle->{{AbsoluteThickness[1],White},{AbsoluteThickness[1],White}},BaseStyle -> {FontFamily->"Arial", FontSize -> 12, PrivateFontOptions->{"OperatorSubstitution"->False,"WindowsUseTrueTypeNames"->True}},ImageSize->650,PlotLabel->"Dre",PlotRange->All]},{
Plot[{
If[Complete,ToConcentration[Are[t]]/.microsim,Null],
If[OrderedSeq,ToConcentration[Are[t]]/.microOrderedSim,Null],
If[RandomSeq, ToConcentration[Are[t]]/.microRandomSim,Null],
If[PingPong, ToConcentration[Are[t]]/.microPingPongSim,Null],
If[Complete,(Are[t])/.macrosim,Null],
If[OrderedSeq,(Are[t])/.macrosimo,Null],
If[RandomSeq,(Are[t])/.macrosimr,Null],
If[PingPong,(Are[t])/.macrosimp,Null]
},{t,0,time},PlotStyle->{Directive[Green,Thick],Directive[Blue,Thick],Directive[Orange,Thick],Directive[Red,Thick],Directive[Black,Thick,Dashed],Directive[Darker[Blue],Thick,Dotted],Directive[Darker[Orange],Thick,Dotted],Directive[Darker[Red],Thick,Dotted]},Frame->{{True,False},{True,False}},FrameLabel->{Style["Time (S)",Bold], Style["Concentration (μM)",Bold]},FrameStyle->{{AbsoluteThickness[1],White},{AbsoluteThickness[1],White}},BaseStyle -> {FontFamily->"Arial", FontSize -> 12, PrivateFontOptions->{"OperatorSubstitution"->False,"WindowsUseTrueTypeNames"->True}},ImageSize->650,PlotLabel->"Are",PlotRange->All],
Plot[{
If[Complete,ToConcentration[Aox[t]]/.microsim,Null],
If[OrderedSeq,ToConcentration[Aox[t]]/.microOrderedSim,Null],
If[RandomSeq, ToConcentration[Aox[t]]/.microRandomSim,Null],
If[PingPong, ToConcentration[Aox[t]]/.microPingPongSim,Null],
If[Complete,(Aox[t])/.macrosim,Null],
If[OrderedSeq,(Aox[t])/.macrosimo,Null],
If[RandomSeq,(Aox[t])/.macrosimr,Null],
If[PingPong,(Aox[t])/.macrosimp,Null]
},{t,0,time},PlotStyle->{Directive[Green,Thick],Directive[Blue,Thick],Directive[Orange,Thick],Directive[Red,Thick],Directive[Black,Thick,Dashed],Directive[Darker[Blue],Thick,Dotted],Directive[Darker[Orange],Thick,Dotted],Directive[Darker[Red],Thick,Dotted]},ImageSize->650,PlotLabel->"Aox",PlotRange->All,Frame->{{True,False},{True,False}},FrameLabel->{Style["Time (S)",Bold], Style["Concentration (μM)",Bold]},FrameStyle->{{AbsoluteThickness[1],White},{AbsoluteThickness[1],White}},BaseStyle -> {FontFamily->"Arial", FontSize -> 12, PrivateFontOptions->{"OperatorSubstitution"->False,"WindowsUseTrueTypeNames"->True}}]}}],{time,1,endTime,1},Row[{Control[{Complete,{True,False}}],Control[{OrderedSeq,{True,False}}],Control[{RandomSeq,{True,False}}],Control[{PingPong,{True,False}}]},Spacer[20]]],TextCell["Green Solid = Complete Micro, Black Dashed = Complete Macro, Blue Solid = Ordered Micro, Blue Dashed = Ordered Macro, Orange Solid = Random Micro, Orange Dashed = Random Macro, Red Solid = Micro PingPong, Red Dashed = Macro PingPong"]},WindowSize->All,WindowTitle->"Comparing Micro/Macro Output-Free Reactants Concentrations Time-Series"]]

Total Substrate Concentration

Concentration time-series composed of both free reactant concentrations (e.g. Dre, Dox, Are, Aox) and reactants bound to the enzyme complex (Dre⋄Epp, etc.).

totalSubstrateCompareButton=Button["Click to Show Comparison of Micro/Macro Output-Total Reactant Concentrations Time-Series",CreateDocument[{Manipulate[Grid[{{
Plot[{
If[Complete,ToConcentration[TotalSubstrate[microsim,Dox,t]],Null],
If[OrderedSeq,ToConcentration[TotalSubstrate[microOrderedSim,Dox,t]],Null],
If[RandomSeq, ToConcentration[TotalSubstrate[microRandomSim,Dox,t]],Null],
If[PingPong, ToConcentration[TotalSubstrate[microPingPongSim,Dox,t]],Null],
If[Complete,(Dox[t])/.macrosim,Null],
If[OrderedSeq,(Dox[t])/.macrosimo,Null],
If[RandomSeq,(Dox[t])/.macrosimr,Null],
If[PingPong,(Dox[t])/.macrosimp,Null]
},{t,0,time},PlotStyle->{Directive[Green,Thick],Directive[Blue,Thick],Directive[Orange,Thick],Directive[Red,Thick],Directive[Black,Thick,Dashed],Directive[Darker[Blue],Thick,Dotted],Directive[Darker[Orange],Thick,Dotted],Directive[Darker[Red],Thick,Dotted]},AxesLabel->{"Time (sec)","Concentration (μM)"},ImageSize->650,PlotLabel->"Dox",PlotRange->All],
Plot[{
If[Complete,ToConcentration[TotalSubstrate[microsim,Dre,t]],Null],
If[OrderedSeq,ToConcentration[TotalSubstrate[microOrderedSim,Dre,t]],Null],
If[RandomSeq, ToConcentration[TotalSubstrate[microRandomSim,Dre,t]],Null],
If[PingPong, ToConcentration[TotalSubstrate[microPingPongSim,Dre,t]],Null],
If[Complete,(Dre[t])/.macrosim,Null],
If[OrderedSeq,(Dre[t])/.macrosimo,Null],
If[RandomSeq,(Dre[t])/.macrosimr,Null],
If[PingPong,(Dre[t])/.macrosimp,Null]
},{t,0,time},PlotStyle->{Directive[Green,Thick],Directive[Blue,Thick],Directive[Orange,Thick],Directive[Red,Thick],Directive[Black,Thick,Dashed],Directive[Darker[Blue],Thick,Dotted],Directive[Darker[Orange],Thick,Dotted],Directive[Darker[Red],Thick,Dotted]},AxesLabel->{"Time (sec)","Concentration (μM)"},ImageSize->650,PlotLabel->"Dre",PlotRange->All]},{
Plot[{
If[Complete,ToConcentration[TotalSubstrate[microsim,Are,t]],Null],
If[OrderedSeq,ToConcentration[TotalSubstrate[microOrderedSim,Are,t]],Null],
If[RandomSeq, ToConcentration[TotalSubstrate[microRandomSim,Are,t]],Null],
If[PingPong, ToConcentration[TotalSubstrate[microPingPongSim,Are,t]],Null],
If[Complete,(Are[t])/.macrosim,Null],
If[OrderedSeq,(Are[t])/.macrosimo,Null],
If[RandomSeq,(Are[t])/.macrosimr,Null],
If[PingPong,(Are[t])/.macrosimp,Null]
},{t,0,time},PlotStyle->{Directive[Green,Thick],Directive[Blue,Thick],Directive[Orange,Thick],Directive[Red,Thick],Directive[Black,Thick,Dashed],Directive[Darker[Blue],Thick,Dotted],Directive[Darker[Orange],Thick,Dotted],Directive[Darker[Red],Thick,Dotted]},AxesLabel->{"Time (sec)","Concentration (μM)"},ImageSize->650,PlotLabel->"Are",PlotRange->All],
Plot[{
If[Complete,ToConcentration[TotalSubstrate[microsim,Aox,t]],Null],
If[OrderedSeq,ToConcentration[TotalSubstrate[microOrderedSim,Aox,t]],Null],
If[RandomSeq, ToConcentration[TotalSubstrate[microRandomSim,Aox,t]],Null],
If[PingPong, ToConcentration[TotalSubstrate[microPingPongSim,Aox,t]],Null],
If[Complete,(Aox[t])/.macrosim,Null],
If[OrderedSeq,(Aox[t])/.macrosimo,Null],
If[RandomSeq,(Aox[t])/.macrosimr,Null],
If[PingPong,(Aox[t])/.macrosimp,Null]
},{t,0,time},PlotStyle->{Directive[Green,Thick],Directive[Blue,Thick],Directive[Orange,Thick],Directive[Red,Thick],Directive[Black,Thick,Dashed],Directive[Darker[Blue],Thick,Dotted],Directive[Darker[Orange],Thick,Dotted],Directive[Darker[Red],Thick,Dotted]},AxesLabel->{"Time (sec)","Concentration (μM)"},ImageSize->650,PlotLabel->"Aox",PlotRange->All]}}],{time,1,endTime,1},Row[{Control[{Complete,{True,False}}],Control[{OrderedSeq,{True,False}}],Control[{RandomSeq,{True,False}}],Control[{PingPong,{True,False}}]},Spacer[20]]],TextCell["Green Solid = Complete Micro, Black Dashed = Complete Macro, Blue Solid = Ordered Micro, Blue Dashed = Ordered Macro, Orange Solid = Random Micro, Orange Dashed = Random Macro, Red Solid = Micro PingPong, Red Dashed = Macro PingPong"]},WindowSize->All,WindowTitle->"Comparing Micro/Macro Output-Total Reactant Concentrations Time-Series"]]

Compare Flux/Rate of Change Time Series

fluxCompareButton=Button["Click to Show Comparison of Flux (Velocity) Time-Series",CreateDocument[{
Manipulate[
Plot[{
If[Complete,(ToConcentration[TotalSubstrateD[microsim,Dre,t]]),Null],
If[Complete,(ToConcentration[TotalSubstrateD[microsim,Aox,t]]),Null],
If[PingPong,(ToConcentration[TotalSubstrateD[microPingPongSim,Dre,t]]),Null],
If[PingPong,(ToConcentration[TotalSubstrateD[microPingPongSim,Aox,t]]),Null],
If[OrderedSeq,(ToConcentration[TotalSubstrateD[microOrderedSim,Dre,t]]),Null],
If[OrderedSeq,(ToConcentration[TotalSubstrateD[microOrderedSim,Aox,t]]),Null],
If[RandomSeq,(ToConcentration[TotalSubstrateD[microRandomSim,Dre,t]]),Null],
If[RandomSeq,(ToConcentration[TotalSubstrateD[microRandomSim,Aox,t]]),Null],
(Dre'[t]/.macrosim)
},{t,0.3,time},
PlotStyle->{Directive[Green,Thick],Directive[Green,Thick,Dashed],Directive[Red,Thick],Directive[Red,Dashed],Directive[Blue,Thick],Directive[Blue,Dashed],Directive[Orange,Thick],Directive[Orange,Dashed],Directive[Black,Thick,Dashed]
},Frame->{{True,False},{True,False}},FrameLabel->{Style["Time (S)",Bold], Style["Flux (((μM)/s))",Bold]},FrameStyle->{{AbsoluteThickness[1],White},{AbsoluteThickness[1],White}},BaseStyle -> {FontFamily->"Arial", FontSize -> 12, PrivateFontOptions->{"OperatorSubstitution"->False,"WindowsUseTrueTypeNames"->True}},ImageSize->700,PlotLabel->"Dre' Aox'",PlotRange->All],{time,1,endTime,1},Row[{Control[{Complete,{True,False}}],Control[{PingPong,{True,False}}],Control[{OrderedSeq,{True,False}}],Control[{RandomSeq,{True,False}}]},Spacer[20]]],TextCell["Solid Color Lines = [Dre], Dashed Color Lines = [Aox], Dashed Black Line = Macro Model"]},WindowSize->All,WindowTitle->"Comparison of Flux (Velocity) Time-Series"]]

Stochastic Simulation Algorithm (xSSA)

Gillespi' s Stochastic Simulation Algorithm is built into another package of xCellerator, xSSA, which can take the same reaction system input format as that of xCellerator. However, because SSA is an exact stochastic simulation, it is very computationally expansive and may quickly become intractable in terms of simulation time and memory storage to record all the time points. Thus, for demonstration purpose, the reaction concentrations, simulation end time, and the order of binding reactions have all been reduced to facilitate shorter simulation time. Note, however, that despite the reduced system, the average of the SSA trials match closely to that of the steady-state chemiosmotic rate equation output.

Micro - Ordered Sequential Bi Bi

endTime2=60; (* Reduced simulated time range - seconds *)
trials=10; (* Number of trials of SSA simulation to average *)

Set2={
complexNumber->ToMolecules[12.5],
DreNumber->ToMolecules[100],
DoxNumber->0,
AoxNumber->ToMolecules[100],
AreNumber->0
};(* Reduced reactant and enzyme concentrations *)
ComplexIb={
kCat->0.3,(*tri-molecular catalysis*)
kCatb->0.3,(*bi-molecular catalysis*)
kCatu->0.3,(*uni-molecular catalysis*)
boundaryTransfer->0.3 (10^3),(* reaction order of the boundary electron transfer reactions *)
BindingOrder->0.3 (10^6)(* The order of the binding steps (M^(-1))*),
Ksd->12 (10^(-6))(*M*),(* Binding constant for donor *)
Ksad->12 (10^(-6))(*M*), (* Binding constant for donor after binding of acceptor *)
Ksa->16.5 (10^(-6))(*M*), (* Binding constant for acceptor *)
Ksda->16.5 (10^(-6))(*M*), (* Binding constant for acceptor after binding of donor *)
Krd->0.01,(* Tri-molecular boundary electron transfer with the donor EC *)
Krdb->0.01,(* Bi-molecular boundary electron transfer with the donor EC *)
Kra->0.1,(* Tri-molecular boundary electron transfer with the acceptor EC *)
Krab->0.1,(* Bi-molecular boundary electron transfer with the acceptor EC *)
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*)
νdr-> 1,
νdo-> 1,
νar-> 1,
νao-> 1
}; (* Reduced binding order *)
microInitialConditions/.Set2

macrosim2= NDSolve[Join[myODE/.adjust/.ComplexIb/.Set2,initial/.ComplexIb/.Set2],var,{t,startTime,endTime2}]; (*Steady-state output based on the reduced parameter values *)

completeSSA := SSA[EC1[(Overscript[{Dre,Dox}⇄{Are,Aox}, E]),ComplexIb,ReactionType->"Complete"], endTime2, microInitialConditions/.Set2, {}]; (* SSA simulation notation for one trial *)

SSAButton=Button["Click to Compare SSA and Steady-State Rate Equation Concentration Time-Series (WARNING! Might take a long time!)",
doc=CreateDocument[{
ExpressionCell[Defer[microSSASims=Table[completeSSA,{trials}];],"Input"],
ExpressionCell[Defer[Cstats=SSAMean[microSSASims];],"Input"],ExpressionCell[Defer[Plot[{ToConcentration[TotalSubstrate[{"Means"/.Cstats},Dre,t]],ToConcentration[TotalSubstrate[{"MeanPlusSigma"/.Cstats},Dre,t]],ToConcentration[TotalSubstrate[{"MeanMinusSigma"/.Cstats},Dre,t]],
ToConcentration[TotalSubstrate[{"MeanPlus2Sigma"/.Cstats},Dre,t]],ToConcentration[TotalSubstrate[{"MeanMinus2Sigma"/.Cstats},Dre,t]],Dre[t]/.macrosim2},{t,0,endTime2} ,PlotStyle->{Directive[Green,Thick,Dashed],Directive[Green,Dashed],Directive[Green,Dashed],Directive[Green,Dotted],Directive[Green,Dotted],Directive[Green,Thick]}, Filling->{1->{{2},Directive[Opacity[0.9],LightGreen]},1->{{3},Directive[Opacity[0.9],LightGreen]},2->{{4},Directive[Opacity[0.3],LightGreen]},3->{{5},Directive[Opacity[0.3],LightGreen]}},PlotRange->{{0,60},{0,100}},ImageSize->800,Frame->{{True,False},{True,False}},FrameLabel->{Style["Time (S)",Bold], Style["Concentration (μM)",Bold]},FrameStyle->{{AbsoluteThickness[1],White},{AbsoluteThickness[1],White}},BaseStyle -> {FontFamily->"Arial", FontSize -> 12, PrivateFontOptions->{"OperatorSubstitution"->False,"WindowsUseTrueTypeNames"->True}}]],"Input"]},WindowTitle->"SSA and Steady-State Rate Equation Concentration Time-Series Comparison"];
SelectionMove[doc,All,Notebook];
SelectionEvaluate[doc];]

Reaction Pathway Composition

Mitochondria Oxidative Phosphorlyation

complexI=EC1[(Overscript[{NADH,NAD}⇄{CoQH2,CoQ}, CI]),ComplexI,ReactionType->"PingPong"];
complexIII=EC1[(Overscript[{CoQH2,CoQ}⇄{cytCred,cytCox}, CIII]),ComplexIII,ReactionType->"Ordered"];
complexIV=EC1[(Overscript[{cytCred,cytCox}⇄{H2O,O2}, CIV]),ComplexIV,ReactionType->"Random"];

reactionSystem=Join[complexI,complexIII,complexIV];
reactionSystem//TableForm

complex1IC={CIpp->3,NADH->ToMolecules[100],NAD->0,CoQ->ToMolecules[90],CoQH2->ToMolecules[10]};
complex3IC={CIIIpp->6,CoQ->ToMolecules[90],CoQH2->ToMolecules[10],cytCox->ToMolecules[25],cytCred->ToMolecules[25]};
complex4IC={CIVpp->9,cytCox->ToMolecules[25],cytCred->ToMolecules[25],O2->ToMolecules[1000],H2O->ToMolecules[1000]};
pathwayIC=Union[complex1IC,complex3IC,complex4IC]

pathwaySim=run[reactionSystem,{0,endTime},rates->{},initialConditions->pathwayIC,frozen->{H2O,O2},plot->False]

freeReactants={NADH,NAD,CoQ,CoQH2,cytCox,cytCred,O2,H2O};
complexISpecies={CIpp,CInp,CIpn,CInn};
complexIIISpecies={CIIIpp,CIIInp,CIIIpn,CIIInn};
complexIVSpecies={CIVpp,CIVnp,CIVpn,CIVnn};
gridPlot[pathwaySim,freeReactants,{0,600}]
gridPlot[pathwaySim,complexISpecies,{0,600}]
gridPlot[pathwaySim,complexIIISpecies,{0,600}]
gridPlot[pathwaySim,complexIVSpecies,{0,600}]

Appendix A - Enumeration of Elementary Reaction Steps Within RedoxMech

Molecular Species/State Variables

Donor={Dre, Dox};
Acceptor={Aox,Are};
Enz={Epp,Enp,Epn,Enn};

Binding Reactions

Donor Electron Carrier Binds First to Enzyme Complex

In the simpliest model, we make the simplification that the two Donor redox species bind equally to the four possible forms of the Enz Complex with a certain equilibrium dissociation constant Ks,d=(kr/kf)=(([D][Enz])/([D⋄Enz])), thus making reactions 1-8 invarient in rate

donorFirst=TwoReactantsReaction[Donor,Enz,D1,ReverseIndex->False];
donorFirstTable=ReactionTable[donorFirst]

Acceptor Electron Carrier Binds First to Enzyme Complex

In the simpliest model, we make the simplification that the two Acceptor redox species bind equally to the four possible forms of the Enz Complex with a certain equilibrium dissociation constant Ks,a=(kr/kf)=(([A][Enz])/([Enz⋄A])), thus making reactions 9-16 invarient in rate

acceptorFirst=TwoReactantsReaction[Enz,Acceptor,A1,ReverseIndex->True];
acceptorFirstTable=ReactionTable[acceptorFirst]

Acceptor Electron Carrier Binds Second to Enzyme Complex

The 16 varients of the ternary D⋄Enz⋄A molecule can be obtained through either sequence of binding of donor or acceptor. The order of the sequence of binding events separates the different mechanisms such as Ordered Sequencial, Random Sequencial, or Ping-Pong. In the simpliest model, we assume parallel binding of the donor and acceptor ECs and thus the secondary binding obeys the same equilibrium dissociation constant as the first.

acceptorSecond=TwoReactantsReaction[donorFirst[[2]],Acceptor,A2,ReverseIndex->True];
acceptorSecondTable=ReactionTable[acceptorSecond]

Donor Electron Carrier Binds Second to Enzyme Complex

The 16 varients of the ternary D⋄Enz⋄A molecule can be obtained through either sequence of binding of donor or acceptor. The order of the sequence of binding events separates the different mechanisms such as Ordered, Random, or Ping-Pong. In the simpliest model, we assume parallel binding of the donor and acceptor ECs and thus the secondary binding obeys the same equilibrium dissociation constant as the first.

donorSecond=TwoReactantsReaction[Donor,acceptorFirst[[2]],D2,ReverseIndex->False];
donorSecondTable=ReactionTable[donorSecond]

Electron Transfer Reactions (Tri-Molecular)

Electron Transfer From Donor EC to the Complex Enzyme Via Rh Reaction Center #49-52

donorToEnzTri={
{{Dre⋄Epp⋄Aox->Dox⋄Enp⋄Aox,(Subscript[kf, TD])[1]},{Dox⋄Enp⋄Aox->Dre⋄Epp⋄Aox,(Subscript[kr, TD])[1]}},
{{Dre⋄Epn⋄Aox->Dox⋄Enn⋄Aox,(Subscript[kf, TD])[2]},{Dox⋄Enn⋄Aox->Dre⋄Epn⋄Aox,(Subscript[kr, TD])[2]}},
{{Dre⋄Epp⋄Are->Dox⋄Enp⋄Are,(Subscript[kf, TD])[3]},{Dox⋄Enp⋄Are->Dre⋄Epp⋄Are,(Subscript[kr, TD])[3]}},
{{Dre⋄Epn⋄Are->Dox⋄Enn⋄Are,(Subscript[kf, TD])[4]},{Dox⋄Enn⋄Are->Dre⋄Epn⋄Are,(Subscript[kr, TD])[4]}}};
donorToEnzTriTable=TransferTable[donorToEnzTri]

Electron Transfer Between the Reaction Centers Rh and Rl of the Complex Enzyme #53-56

EnzRhToRlTri=
{{{Dre⋄Enp⋄Aox->Dre⋄Epn⋄Aox,(Subscript[kf, TT])[1]},{Dre⋄Epn⋄Aox->Dre⋄Enp⋄Aox,(Subscript[kr, TT])[1]}},
{{Dox⋄Enp⋄Aox->Dox⋄Epn⋄Aox,(Subscript[kf, TT])[2]},{Dox⋄Epn⋄Aox->Dox⋄Enp⋄Aox,(Subscript[kr, TT])[2]}},
{{Dre⋄Enp⋄Are->Dre⋄Epn⋄Are,(Subscript[kf, TT])[3]},{Dre⋄Epn⋄Are->Dre⋄Enp⋄Are,(Subscript[kr, TT])[3]}},
{{Dox⋄Enp⋄Are->Dox⋄Epn⋄Are,(Subscript[kf, TT])[4]},{Dox⋄Epn⋄Are->Dox⋄Enp⋄Are,(Subscript[kr, TT])[4]}}};
EnzRhToRlTriTable=TransferTable[EnzRhToRlTri]

Electron Transfer From the Complex Enzyme to the Acceptor EC Via Rl Reaction Center #57-60

EnzToAcceptorTri=
{{{Dre⋄Enn⋄Aox->Dre⋄Enp⋄Are,(Subscript[kf, TA])[1]},{Dre⋄Enp⋄Are->Dre⋄Enn⋄Aox,(Subscript[kr, TA])[1]}},
{{Dox⋄Enn⋄Aox->Dox⋄Enp⋄Are,(Subscript[kf, TA])[2]},{Dox⋄Enp⋄Are->Dox⋄Enn⋄Aox,(Subscript[kr, TA])[2]}},
{{Dre⋄Epn⋄Aox->Dre⋄Epp⋄Are,(Subscript[kf, TA])[3]},{Dre⋄Epp⋄Are->Dre⋄Epn⋄Aox,(Subscript[kr, TA])[3]}},
{{Dox⋄Epn⋄Aox->Dox⋄Epp⋄Are,(Subscript[kf, TA])[4]},{Dox⋄Epp⋄Are->Dox⋄Epn⋄Aox,(Subscript[kr, TA])[4]}}};
EnzToAcceptorTriTable=TransferTable[EnzToAcceptorTri]

Electron Transfer Reactions (Bi-Molecular)

Electron Transfer From Donor EC to the Complex Enzyme Via Rh Reaction Center

donorToEnzBi=
{{{Dre⋄Epp->Dox⋄Enp,(Subscript[kf, BD])[1]},{Dox⋄Enp->Dre⋄Epp,(Subscript[kr, BD])[1]}},
{{Dre⋄Epn->Dox⋄Enn,(Subscript[kf, BD])[2]},{Dox⋄Enn->Dre⋄Epn,(Subscript[kr, BD])[2]}}};
donorToEnzBiTable=TransferTable[donorToEnzBi]

Electron Transfer Between the Reaction Centers Rh and Rl of the Complex Enzyme

EnzRhToRlBi=
{{{Dre⋄Enp->Dre⋄Epn,(Subscript[kf, BT])[1]},{Dre⋄Epn->Dre⋄Enp,(Subscript[kr, BT])[1]}},
{{Dox⋄Enp->Dox⋄Epn,(Subscript[kf, BT])[2]},{Dox⋄Epn->Dox⋄Enp,(Subscript[kr, BT])[2]}},
{{Enp⋄Aox->Epn⋄Aox,(Subscript[kf, BT])[3]},{Epn⋄Aox->Enp⋄Aox,(Subscript[kr, BT])[3]}},
{{Enp⋄Are->Epn⋄Are,(Subscript[kf, BT])[4]},{Epn⋄Are->Enp⋄Are,(Subscript[kr, BT])[4]}}
};
EnzRhToRlBiTable=TransferTable[EnzRhToRlBi]

Electron Transfer From the Complex Enzyme to the Acceptor EC Via Rl Reaction Center

EnzToAcceptorBi=
{{{Epn⋄Aox->Epp⋄Are,(Subscript[kf, BA])[1]},{Epp⋄Are->Epn⋄Aox,(Subscript[kr, BA])[1]}},
{{Enn⋄Aox->Enp⋄Are,(Subscript[kf, BA])[2]},{Enp⋄Are->Enn⋄Aox,(Subscript[kr, BA])[2]}}};
EnzToAcceptorBiTable=TransferTable[EnzToAcceptorBi]

Electron Transfer Reactions (Uni-Molecular)

Electron Transfer Between the Reaction Centers Rh and Rl of the Complex Enzyme

EnzRhToRlUni=
{{{Enp->Epn,(Subscript[kf, UT])[1]},{Epn->Enp,(Subscript[kr, UT])[1]}}
};
EnzRhToRlUniTable=TransferTable[EnzRhToRlUni]

Reaction Groups

BindingReactions=Flatten[{donorFirst[[1]],acceptorFirst[[1]],acceptorSecond[[1]],donorSecond[[1]]},2];

BindingHeadings=Join[
Flatten[Table[Flatten[{"D First (molecule/s)","  Dissociation (1/s)"}],{i,1,8}]],
Flatten[Table[Flatten[{"A First (molecule/s)","  Dissociation (1/s)"}],{i,1,8}]],
Flatten[Table[Flatten[{"A Second(molecule/s)","  Dissociation (1/s)"}],{i,1,16}]],
Flatten[Table[Flatten[{"D Second(molecule/s)","  Dissociation (1/s)"}],{i,1,16}]]
];

TernaryTransfers=Flatten[{donorToEnzTri,EnzRhToRlTri,EnzToAcceptorTri},2];

BinaryTransfers=Flatten[{donorToEnzBi,EnzRhToRlBi,EnzToAcceptorBi},2];

UninaryTransfers=Flatten[{EnzRhToRlUni},2];

Reaction Mechanism

(*microReactions=Join[BindingReactions,TernaryTransfers,BinaryTransfers,UninaryTransfers];*)

(*OrderedSeqReactions=Join[
PickReactions[Flatten[donorFirst[[1]],1],{1}],PickReactions[Flatten[acceptorFirst[[1]],1],{5}],PickReactions[Flatten[acceptorSecond[[1]],1],{1}],PickReactions[Flatten[donorSecond[[1]],1],{13}],PickReactions[TernaryTransfers,{1,6,12}],{{Enn->∅,0},{Enp->∅,0},{Epn->∅,0}}
]*)

(*RandomSeqReactions=Join[PickReactions[Flatten[donorFirst[[1]],1],{1,5}],PickReactions[Flatten[acceptorFirst[[1]],1],{1,5}],PickReactions[Flatten[acceptorSecond[[1]],1],{1,13}],PickReactions[Flatten[donorSecond[[1]],1],{1,13}],PickReactions[TernaryTransfers,{1,6,12}],{{Enn->∅,0},{Enp->∅,0},{Epn->∅,0}}
]*)

(*PingPongReactions=Join[
PickReactions[Flatten[donorFirst[[1]],1],{1,6}],
PickReactions[Flatten[acceptorFirst[[1]],1],{2,5}],
PickReactions[BinaryTransfers,{1,5,7}],{{Enn->∅,0},{Epn->∅,0}}
]*)

Parameter Specification - Complete Dynamics

donorFirstParam=Flatten[Table[{
(Subscript[kf, D1])[i]->N[(BindingOrder/(nA Volume))],
(Subscript[kr, D1])[i]->N[BindingOrder Ksd]},{i,1,8}]];
acceptorFirstParam=Flatten[Table[{
(Subscript[kf, A1])[i]->N[(BindingOrder/(nA Volume))],
(Subscript[kr, A1])[i]->N[BindingOrder Ksa]},{i,1,8}]];
acceptorSecondParam=Flatten[Table[{
(Subscript[kf, A2])[i]->N[(BindingOrder/(nA Volume))],
(Subscript[kr, A2])[i]->N[BindingOrder Ksda]},{i,1,16}]];
donorSecondParam=Flatten[Table[{
(Subscript[kf, D2])[i]->N[(BindingOrder/(nA Volume))],
(Subscript[kr, D2])[i]->N[BindingOrder Ksad]},{i,1,16}]];
donorToEnzTriParam=Flatten[Table[{
(Subscript[kf, TD])[i]->boundaryTransfer,
(Subscript[kr, TD])[i]->boundaryTransfer Krd},{i,1,4}]];
EnzRhToRlTriParam=Flatten[Table[{
(Subscript[kf, TT])[i]->kCat,
(Subscript[kr, TT])[i]->N[(kCat/Keq)]},{i,1,4}]];
EnzToAcceptorTriParam=Flatten[Table[{
(Subscript[kf, TA])[i]->boundaryTransfer,
(Subscript[kr, TA])[i]->boundaryTransfer Kra},{i,1,4}]];
donorToEnzBiParam=Flatten[Table[{
(Subscript[kf, BD])[i]->boundaryTransfer,
(Subscript[kr, BD])[i]->boundaryTransfer Krdb},{i,1,2}]];
EnzRhToRlBiParam=Flatten[Table[{
(Subscript[kf, BT])[i]->kCatb,
(Subscript[kr, BT])[i]->N[(kCatb/Keq)]},{i,1,4}]];
EnzToAcceptorBiParam=Flatten[Table[{
(Subscript[kf, BA])[i]->boundaryTransfer,
(Subscript[kr, BA])[i]->boundaryTransfer Krab},{i,1,2}]];
EnzRhToRlUniParam=Flatten[Table[{
(Subscript[kf, UT])[i]->kCatu,
(Subscript[kr, UT])[i]->N[(kCatu/Keq)]},{i,1,1}]];
completeParam = Flatten[{donorFirstParam, acceptorFirstParam, acceptorSecondParam, donorSecondParam, donorToEnzTriParam, EnzRhToRlTriParam, EnzToAcceptorTriParam, donorToEnzBiParam, EnzRhToRlBiParam, EnzToAcceptorBiParam, EnzRhToRlUniParam}];
reactionListButton=Button["Click to Show Complete Reaction/Parameter List",CreateDocument[{TextCell["Complete Reaction Set","Section"],
TableForm[lowLevelReactions[BindingReactions],TableHeadings->{BindingHeadings,{"Binding/Dissociation Reaction","Rate"}}],
TableForm[lowLevelReactions[TernaryTransfers],TableHeadings->{Flatten[Table[{"Forward","Reverse"},{i,1,12}]],{"Ternary Complex Electron Transfer","Rate (1/s)"}}],
TableForm[lowLevelReactions[BinaryTransfers],TableHeadings->{Flatten[Table[{"Forward","Reverse"},{i,1,8}]],{"Binary Complex Electron Transfer","Rate (1/s)"}}],
TableForm[lowLevelReactions[UninaryTransfers],TableHeadings->{{"Forward","Reverse"},{"Uninary Complex Electron Transfer","Rate (1/s)"}}],TextCell["Ordered Sequential Reaction Set","Section"],
TableForm[lowLevelReactions[OrderedSeqReactions],TableHeadings->{Flatten[Table[{"Forward","Reverse"},{i,1,7}]],{"Reaction","Rate"}}],TextCell["Random Sequential Reaction Set","Section"],
TableForm[lowLevelReactions[RandomSeqReactions],TableHeadings->{Flatten[Table[{"Forward","Reverse"},{i,1,11}]],{"Reaction","Rate"}}],
TextCell["Ping-Pong Reaction Set","Section"],
TableForm[lowLevelReactions[PingPongReactions],TableHeadings->{Flatten[Table[{"Forward","Reverse"},{i,1,7}]],{"Reaction","Rate"}}]},WindowTitle->"Reaction/Parameter List"]]

Appendix B - RedoxMech Feature Contents

The following code automatically creates a separate dialog window containing all the features of the RedoxMech package

CreateDialog[{TextCell["Available Dynamic Contents"],reactionListButton,pingPongReportButton,orderedReportButton,randomReportButton,freeSubstrateCompareButton,totalSubstrateCompareButton,fluxCompareButton,SSAButton,TextCell["Change the number of trials for SSA (before starting the simulation)"],InputField[Dynamic[trials],Number],DefaultButton[]},WindowTitle->"RedoxMech Features"];