A Voltammetric Perspective of Multi-Electron and Proton Transfer in Protein Redox Chemistry: Insights From Computational Analysis of Escherichia coli HypD Fourier Transformed Alternating Current Voltammetry

This paper explores the impact of pH on the mechanism of reversible disulfide bond (CysS-SCys) reductive breaking and oxidative formation in Escherichia coli hydrogenase maturation factor HypD, a protein which forms a highly stable adsorbed film on a graphite electrode. To achieve this, low frequency (8.96 Hz) Fourier transformed alternating current voltammetric (FTACV) experimental data was used in combination with modelling approaches based on Butler-Volmer theory with a dual polynomial capacitance model, utilizing an automated two-step fitting process conducted within a Bayesian framework. We previously showed that at pH 6.0 the protein data is best modelled by a redox reaction of two separate, stepwise one-electron, one-proton transfers with slightly “crossed” apparent reduction potentials that incorporate electron and proton transfer terms (Eapp20 > Eapp10). Remarkably, rather than collapsing to a concerted two-electron redox reaction at more extreme pH, the same two-stepwise one-electron transfer model with Eapp20 > Eapp10 continues to provide the best fit to FTACV data measured across a proton concentration range from pH 4.0 to pH 9.0. A similar, small level of crossover in reversible potentials is also displayed in overall two-electron transitions in other proteins and enzymes, and this provides access to a small but finite amount of the one electron reduced intermediate state.


Supplementary Material for "A Voltammetric Perspective of Multi-Electron and Proton Transfer in Protein Redox Chemistry: Insights from Computational Analysis of Escherichia coli HypD Fourier Transformed Alternating Current Voltammetry"
Contents 1 H2ase-MFHypD DC and FTACV data at pH 6.0 S2 2 Nernst plot S3 3 H2ase-MFHypD FTACV data from pH 4.0 to pH 9.0 S4 4 Mathematical model S5 4.1 Model for two-stepwise one-electron transfers 4.1.1 Dimensionless equations 4.2 Model for one-concerted two-electron transfer 5 Computational methods S11 5. Comparing best fit objective function values for the two models S18 7 Comparing pH 5.0, 6.0, 7.0 and 8.0 experimental data and best fit plots from both reaction models S19 8 Best fit parameter values for the one-concerted two-electron transfer model S20 9 Best fit capacitance parameter values S21 10 Misspecification of the n3 = 2 model S23 11 Verification of H2ase-MFHypD Signal S25 12 Inference and MCMC for Eapp versus pH S26 13 Insensitivity of Faradaic fitting to Uncompensated Resistance S31 14 References S32 S2 1 H2ase-MFHypD DC and FTACV data at pH 6.0

Figure S1
Comparison of pH 6.0 H2ase-MFHypD (a) DCV at a scan rate of 89 mV s -1 and (b) FTACV data obtained with an underlying DC scan rate of 22.4 mV s -1 and a sinewave perturbation of amplitude 150 mV and frequency of 8.96 Hz, where the potential-time methodologies are contrasted as well as the resulting current-potential output voltammogram of DCV and 6 th -10 th harmonics of FTACV (grey letters indicate harmonic number). The vertical red dotted line indicates how the reversible midpoint potential associated with the overall redox reaction (Emid) is extracted from each experiment. Nernst plot Figure S2 Plots derived from the Nernst equation showing how under equilibrium conditions the amount of an "intermediate" state generated during two-stepwise, reversible one-electron transfers depends on the difference in the reversible potentials governing each one-electron process. "Intermediate" is formed by adding one-electron (n1 = 1) to the "oxidised" state (equilibrium potential = 1 0 ) and it is consumed by a one-electron reduction (n2 = 1) that produces the "reduced" state (equilibrium potential = 2 0 ). For comparison, plots for a single one-electron "n = 1" reversible reduction process (equilibrium potential = 1 0 ) and simultaneous "n3 = 2" process (equilibrium potential = 3 0 ) are also included, the latter is indistinguishable from the two-stepwise one-electron transfers when 1 0 is sufficiently more negative than 2 0 .

Figure S3
Stack-plot comparison of 6 th harmonic H2ase-MFHypD FTACV experimental data collected previously (Adamson et al., 2017b). FTACV experiments were conducted with dc scan rate (ν) = 22.4 mV s -1 , AC sinewave frequency (f) = 8.96 Hz, and AC sinewave amplitude ( ) = 150 mV. Experimental pH values are as indicated by the legend. (b) Scatter plot of the midpoint potentials extracted from experiments shown in (a) and reference lines indicating (blue) -59 mV pH -1 gradient extrapolated from the pH 6.0 experimental data point and (yellow) -29 mV pH -1 gradient extrapolated from the pH 9.0 experimental data point.

Mathematical model
The mathematical model used is derived from that described in our previous work on H2ase-MFHypD (Adamson et al., 2017b), with changes to the capacitive current calculations (vide infra). As described in the main paper, two different reaction regimes are considered, both a two stepwise one-electron transfer reaction (n1=n2=1), where a surface-confined intermediate oxidation state is formed, Y (surf) , as in Reaction S1 and Reaction S2 (Reactions 1 and 2 in main paper); and a simultaneous two-electron transfer reaction (n3=2) to convert the surfaceconfined oxidised species (X (surf) ) into a surface-confined reduced species (Z (surf) ), as in Reaction S3 (Reaction 3 in the main paper). In all cases, i red is the potential-dependent reductive/forward rate constant for reaction i, and i ox is the potential-dependent oxidative/backwards rate constant for reaction i.
The "app" subscript is used to denote the fact that conditions cannot be achieved where protonation is absent, so that true E 0 values for non-proton coupled electron transfer processes are unknown, as are protonation equilibrium Ka values. Instead, it is assumed that protonation reactions are reversible (diffusion controlled), and modelling is therefore undertaken by combining the E o , Ka and pH terms into an apparent "E o app" value which is defined as the reversible potential (Bond, 2002;Adamson et al., 2017b). This is possible because the processes maintain the same shape and are simply shifted to less positive potentials by the reversible chemical reactions following the electron transfer step. On this basis, the heterogeneous charge transfer rate constant k 0 app and charge transfer coefficient α 0 app are the electrode kinetic parameters relevant to E 0 app.
Using Butler Volmer theory (Bard and Faulkner, 2000), the forwards i red ( ) and backwards i ox ( ) potential dependent reaction rate constants for an electron transfer step at time are written as: where the subscript i denotes the reaction and the superscript 0 indicates stated conditions, app i 0 is the electron transfer rate coefficient, i is the charge transfer coefficient, is the Faraday constant, is the gas constant, is temperature in Kelvin, i is the number of S6 electrons transferred in the reaction, and r ( ) is the real (effective) potential seen by the working electrode at time . r ( ), differs from the input applied potential ( ) due to uncompensated resistance u resulting in an 'ohmic drop' ( u ). In electrochemistry ohmic drop is the result of the potential induced by the resistance of the electrolyte or any other interface such as surface films or connectors. r ( ) can therefore be formulated as: where tot ( ) is the total current measured and is described later.
In FTACV the input potential as a function of time, ( ), is formulated as Equation S4, the same equation is also given in the main paper: The first term of Equation S4, dc ( ), describes dc contributions to the overall input potential, while the second term, Δ sin( + ), describes the ac contribution; Δ is the amplitude of the sine wave, is the angular frequency ( = 2 , where is frequency in Hz), and is the phase. Thus, DCV can be accommodated by Equation S4, as the second term will be zero.
Furthermore, dc ( ) is formulated accordingly to accommodate cyclic voltammetry: where is the dc scan rate (i.e. rate of change of the linear dc ramp), start is the initial applied potential, and reverse is the potential at which the direction of potential change reverses.

Model for two-stepwise one-electron transfers
To model Reaction S1 and Reaction S2, we let the proportion of surface-confined species in state X be X , the proportion of surface-confined species in state Y be Y , and the proportion of surface-confined species in state Z be Z , with initial conditions X (0) = 1, Y (0) = 0, and Z (0) = 0. The ordinary differential equations governing their behaviour are given by: Note the substitution Y = 1 − X − Z has been used and d Y d has been omitted as , it is therefore not independent and the system may be fully described without it.
There are two contributions to the total current, tot , a capacitive term ( c ) and a Faradaic term ( f ), yielding: where c is defined in units of amps per unit area, as is standard, and is therefore scaled by electrode area (assumed to equal the geometric area, 0.03 cm 2 ). The Faradaic term is due to the electrons transferred between the enzyme and the electrode, and is therefore formulated: where is the surface coverage of enzyme per unit area. The electrode used in the experiment lacks a reliable model for the capacitive term, therefore the following third order polynomial in r is used: where dl , dl1 , dl2 , and dl3 are best fit constants. In this work, different constants dl , dl1 , dl2 , and dl3 are defined for the regions of the experiment either side of the reverse value, this is described by Equation S12 to Equation S15.
The new capacitance model provides a superior fit to the background capacitative current when compared to the model used in our previous work which assumed the background current was independent of scan direction (Adamson et al., 2017b). This is attributed to two reasons: firstly, the new capacitance model incorporates more parameters, therefore it has more degrees of freedom. Secondly, the new model accounts for the physical reality that in our system we would not expect capacitance to be independent of scan direction. If we consider the negative DC current scan direction, ideally complete reduction of the surface confined species will have been achieved by the end of this voltage sweep, resulting in an electrode surface modified by the presence of the reduced species. The presence of a reduced surface-confined species can provide a different contribution to the background current to that encountered with surface-confined species prior to reduction. This potential dependent surface state will lead to hysteresis in the background current, with the profile varying with scan direction. The use of a different models to describe the background when sweeping the potential in the negative and positive directions accommodates this hysteresis.

Dimensionless equations
To aid computation, Equation S7 to Equation S9 were represented in a dimensionless form by substituting the dimensionless variables = tot / 0 , = / 0 , and = / 0 for all potentials and setting: Equation S16 0 is an approximation of electrode coverage with units of moles per unit area. The time dependent variables X and Z are already dimensionless and therefore remain. For Equation S7 this gives: where 1 0 = app1 0 0 and is the dimensionless reaction kinetics parameter. For Equation S8 it gives: where the dimensionless reaction kinetics parameter is 2 0 = app2 0 0 . While for Equation S9 it yields: where the dimensionless electrode coverage = / 0 0 = / 0 ( is introduced to adjust the surface coverage ); and the dimensionless double layer capacitance parameters are:

Equation S23
Additionally: Where the dimensionless uncompensated resistance = u 0 / 0 , r = r / 0 , and any ac input signal is parameterized by its dimensionless amplitude = / 0 , angular frequency = 0 , and phase . Note the conditions on Equation S20 to Equation S26 depend on the sign of 0 ; they are therefore defined above for positive as this is the standard why of defining such equations, however negative and therefore equations with the relative equality changes were used in this work.

Model for one-concerted two-electron transfer
To model Reaction S3, capacitance is as before (Equation S11), and 3 red and 3 ox are calculated from Equation S1 and Equation S2, noting that n = 2. Defining X and Z as the proportion of the X and Z respectively with initial conditions X (0) = 1 and Z (0) = 0 the ordinary differential equation governing the systems behaviour is given by: Note the substitution Z = 1 − X has been used and , it is therefore not independent and the system may be fully described without it.
The Faradaic term is due to the electrons transferred between the enzyme and the electrode, and is therefore formulated as: The dimensionless equations using the same parameters as before (see Dimensionless equations) but with the non-dimensionality constants: are therefore: and: The same two-step fitting approach described in our previous H2ase-MFHypD work was used, 1 first fitting for capacitance in the time domain in regions of the current trace with little or no Faradaic contribution, then fitting the Faradaic parameters in the frequency domain using Fourier transformation. Three models were therefore written, one to simulate and fit only capacitance-current using the equations outlined above and setting f to zero, and two full models to simulate both current components and fit for Faradaic parameters in Fourier space using either a model for two-stepwise one-electron transfers ( 1 = 2 = 1) or oneconcerted two-electron transfer ( 3 = 2). The known experimental parameters are as summarised in Figure S4.

Figure S4
Known parameters of all experimental data, pH in the formulae of start and reverse refer to the pH of the experimental data, 1 and 2 are part of the two sequential one electron transfer model ( 1 = 2 = 1), while 3 is part of the concerted two electron transfer model ( 3 = 2).
Prior to fitting, the raw experimental data used for this work was manipulated in two ways: it was down-sampled so that approximately 200 measurements occurred per period of the applied AC current as is standard practice (~217 was used); and it was ensured that the measurements finished on an integer number of periods as this guarantees the data begins and ends at the same location with respect to the oscillating input signal, easing Fourier transformation. A single period was therefore dropped from the end of all experimental measurements, this had the added benefit of removing larger noisy current spikes present at the end of some of the data. As the same manipulation was done to each experimental current trace fair comparison can be done between experiments. Note that a higher sample rate of ~543 measurements per period was used for plots.

Forward problem
To solve and code these models the dimensionless equations must be solved for the unknown time dependent variables X , Z , and . The same approach was taken as used in our previous work. Firstly, the time derivatives were discretised using the backwards Euler method. These discretised equations were then substituted into Equation S19/Equation S31, giving an implicit equation in terms of X n+1 , X n , Z n+1 , Z n , n+1 , and n . This was then solved using the Newton-Raphson method which was implemented in Python, thus solving the forward problem. An equivalent approach was taken for the 3 = 2 model. Note that the only discretization required for the capacitance model is d d , as only Equation S11 or the first term of Equation S19/Equation S31 when non-dimensionalized is needed for the capacitance only model. S13

Inverse problem
As previously mentioned, a two-step fitting method, very similar to that used in our previous work on H2ase-MFHypD (Adamson et al., 2017b), was used to fit the experimental data. Firstly, background capacitive current was fitted using a capacitance only model and regions of the experimental data which could be safely described using this model. Subsequently, Faradaic current was fitted for using a Fourier transform to separate the Faradaic signal from the data. This process is outlined below. Note the optimisation in this work was carried out in a Bayesian framework, the mathematics of this has been covered in the electrochemistry literature (Gavaghan et al., 2018;Gundry et al., 2021), with extensive further cover available in the statistical literature (Lambert, 2018), and is therefore not repeated in this work.

Finding best fit parameter values for the non-Faradaic capacitance-current model
Fitting of the background capacitance was carried out in the time domain, using the regions of the experimental data that could be described using only the capacitance model. These regions were identified by visual inspection of the current trace and Figure S5 shows them highlighted in blue for a sample of pH 6.0 experimental data. These regions represent the same time-window across all the data because in the Experimental data collection start and rev were adjusted with pH so that the Faradaic region occupied the same time-window at all pH's (see Figure S3). The capacitance parameters vdl , vdl1 , vdl2 , vdl3 , −vdl , −vdl1 , −vdl2 , and −vdl3 were fitted for as well as , , and u . The phase shift, found for capacitance was not used for the Faradaic fitting and was refitted for, this is standard practise due to the sensitivity of the model to the phase.

Figure S5
The current trace for pH 6.0 experiment 3, shaded in blue are the assumed capacitive regions used to fit the capacitance parameters and uncompensated resistance, and shaded in black is the voltage region within which This inverse problem was addressed using the infrastructure provide by the PINTS Python library (Clerx et al., 2019). The model was wrapped as a 'Pints.ForwardModel' to allow use of the automated methods of parameter fitting implemented in the PINTS Python library (Clerx et al., 2019). The problem was then optimised using the CMA-ES optimiser, a uniform prior with box bound parameter given in Figure S6 (transformed using 'Pints.RectangularBoundriesTranformation' to aid fitting), and sum of square difference objective function, ℓ , between simulated capacitance current, c , and experimental current, e , as was done previously (Adamson et al., 2017b).

Equation S32
For each experimental dataset 10 optimisations were performed using the box bound parameters described in Figure S6, the objective function given in Equation S32, and the CMA-ES optimiser from the PINTS Python library (Clerx et al., 2019). A maximum number of iterations of 10,000, or 200 consecutive unchanged iterations with a threshold of 1e-11 being deemed as unchanged, was used to terminate each optimisation. The fit giving the lowest objective functions score was deemed the best and was generated a minimum of 3 times in each experiment. The results of this are presented in Table S4.

Figure S6
Box bound prior distribution of capacitance parameters for inference, with 0 = 2 .

Finding best fit parameter values for the two different reaction models
Fourier transformation is a very natural way of transforming the output of an AC voltammetry experiment, as the constant frequency AC nature of the input signal gives a natural timescale to the output. Indeed it has been shown in the literature that fourth and higher order harmonics form PF-FTACV have negligible capacitive current and therefore provide access to an essentially background free measurement of Faradaic current derived from fast protein redox processes (Lee et al., 2011;Stevenson et al., 2012;Adamson et al., 2017a). As in our previous work on H2ase-MFHypD data we took advantage of this separation using a two-step fitting process (Adamson et al., 2017b), firstly fitting capacitance parameters as outlined above, before determining the Faradaic parameters that define protein redox processes ( i 0 , i 0 , , and ) by fitting to the real and imaginary contributions of the 4th to 12th order harmonics of the ac data using data optimisation methods. Harmonics 4 to 12 are used as harmonics below the 4th describe capacitance and harmonics above the 12 th (for our data) become masked by experimental noise.
Fourier transformation of the experimental and simulated data was therefore carried out using the fast Fourier transform available in the NumPy Python library (Harris et al., 2020). Additional manipulations of the experimental data and simulated data was then carried out to separate out the 4 th to 12 th harmonics. The Fourier transformed data corresponding to 3.5 to 12.5 was used to 'top hat' filter for the harmonics, this generous range either side of the peripheral harmonics was used to account for any phase shift and/or inaccuracy in the experimental reported frequency.
The infrastructure provide by the PINTS Python library was again employed for optimisation, and the model was therefore wrapped as a 'Pints.ForwardModel' (Clerx et al., 2019). The 0≤ dl ≤1e-3 , S15 problem was then optimised using the C-MAES optimiser, a uniform prior with box bound parameter given in Figure S7 (transformed using 'Pints.RectangularBoundriesTranformation' to aid fitting), and an Euclidean distance objective function, ℓ , between simulated Fourier transformed current, n ( k ), and Fourier transformed experimental current, e ( k ): Here k is frequency, and note the use of the modulus operator as n ( k ) and e ( k ) are complex numbers, this transforms the resultant complex number into a real number describing the magnitude of the complex number allowing minimisation of the difference. This differs to the objective function used in our previous work in three ways: in this work a 'top hat' filter of range 3.5 to 12.5 was used to separate harmonics 4-12 and the entire range was used in the objective function as opposed to applying a Kaiser window for each individual harmonic; the objective function is normalised by N; and no windowing functions were used.

Figure S7
Box bound uniform distributions of Faradaic parameters for inference, where min = starts + 0.375( reverse − start ), min = starts + 0.575( reverse − start ), and = / 0 as before. On the left are the 6 parameters fitted for when using the two-sequential one-electron transfer model ( 1 = 2 = 1) and on the right are the 4 parameters fitted for when using the concerted two-electron transfer model ( 3 = 2). An example of the fitting range described by min and max can be seen in Figure S5 (black trace). 10 optimisations were performed using the box bound parameters described in Figure S7, the objective function given in Equation S33, and the CMA-ES optimiser from the PINTS Python library for each experimental current trace (Clerx et al., 2019). A maximum number of iterations of 10,000 or 200 consecutive unchanged iterations with a threshold of 1e-11 being deemed as unchanged was used to terminate each optimisation. The fit giving the lowest objective functions score was deemed the best and was generated a minimum of 4 times in each experiment. The results of this are discussed in the main paper, with results of fitting to the n1=n2=1 model given in Table 1, and results of fitting to the n=2 model given in Table S3. The best fit capacitance parameters previously determined in the first fitting step were used for each data set. Note the phase parameter, , from capacitance fitting was not used and was fitted for again due to a shift that occurs between capacitive and Faradaic S16 regions, as was done in our previous work (Adamson et al., 2017b). The box bound parameters described in Figure S7 were chosen to be consistent with our previous work on H2ase-MFHypD (Adamson et al., 2017b). However, a narrower fitting range for 1 0 , 2 0 , and 3 0 was used as this can be seen to be physically sensible in Figure S5 (black shaded region). S17

Verification of computational methods
All models and inference methods used in this study were verified before use on experimental data. This section details the verification process.

Forward problem verification
The solutions to the 3 = 2 and 1 = 2 = 1 models detailed in Forward Problem were verified by taking an additional approach involving recasting Equation S19/Equation S31 as an ordinary differential equation (ODE) and solving the resulting system of ODE's using the ODE solver 'odeint' from the SciPy Python library using default tolerances -positive error weights of 1. 49E-8 (Virtanen et al., 2020). The two implementations of the models agreedweighted mean square differences of order 1.0E-16 A 2 -across all parameter sets tested verifying that the Newton-Raphson model was correctly implemented. Note the Newton-Raphson approach is preferred for electrochemical problems as it is substantially faster to solve. Additionally, the forward model was verified by reproducing result figures from the literature (Stevenson et al., 2012).

Inverse problem verification
The implementation of the objective function and inference methods for the capacitance and faradaic models were verified by fitting to synthetic data using the experimental relevant parameters given in Figure S4, capacitance and faradaic parameters within the ranges given in Figure S6 and Figure S7 respectively. Both capacitance and faradaic models could reproduce all parameters within the range shown in Figure S6 and Figure S7 respectively without experimental noise. Faradaic parameters were regenerated exactly (to a minimum of 11 S.F) giving an objective function score of ~0.0 (1E-8). Capacitance parameters were fitted for to synthetic data generated using the capacitance model allowing exact recover (to a minimum of 13 S.F) giving an objective function score of ~0.0 (1E-19). Additionally, capacitance parameters were fitted for to synthetic data generated using the full faradaic model to better emulate fitting to experimental data; an objective function of zero was therefore not achieved but the capacitance parameters were regenerated to 2-11 S.F; the fitted capacitance parameters were then used to fit the faradaic parameters yielding an objective function score of 47.8 and parameters accurate to 3-5 S.F. Thus verifying the implementations were working correctly, and preforming as expected. S18

Comparing best fit objective function values for the two models
A lower objective function score (Equation S33) indicates a better fit between a simulated dataset and experimental data. As the two different reaction models (n1=n2=1 versus n3=2) use different non-dimensionalization constants, the objective function scores generated from the "best fit" parameter sets are converted into their dimensional form and then compared. Table S1 shows scores for pH 4.0, pH 5.0, and pH 6.0 data while Table S2 shows the equivalent values for pH 7.0, pH 8.0, and pH 9.0. This analysis emphasises that the n1=n2=1 model provides a better fit to the data for all experiments at all pH's studied.

Table S1
Dimensional objective function values (Equation S33) for the best fit parameter values for the simultaneous two-electron transfer model (n3=2, see Table S3) and the twosequential one-electron transfer model (n1=n2=1 , Table 1) for three ν = 22.4 mV s -1 , f = 8.96 Hz, ΔE = 150 mV H2ase-MFHypD FTACV experiments at pH 4.0, 5.0, and 6.0. The objective function has been reported in dimensional units to allow comparison between the two models. The better fit is shaded in for each experiment.    Table S3) and the twosequential one-electron transfer model (n1=n2=1 , Table 1) for three ν = 22.4 mV s -1 , f = 8.96 Hz, ΔE = 150 mV H2ase-MFHypD FTACV experiments at pH 7.0, 8.0, and 9.0. The objective function has been reported in dimensional units to allow comparison between the two models. The better fit is shaded in for each experiment. Comparing pH 5.0, 6.0, 7.0 and 8.0 experimental data and best fit plots from both reaction models Figure S8 Overlay current-potential plots of (a) pH 5.0, (b) pH 6.0, (c) pH 7.0, and (d) pH 8.0 (grey line) data from experiment 3, and best fit simulated data using (red line) a concerted two-electron (n3 = 2) electron transfer model and (blue line) a stepwise two-sequential oneelectron transfer reaction model (n1=n2=1) for (top) 6 th , (middle) 7 th , and (bottom) 8 th harmonics. The parameter values used to generate the simulations are reported in Table 1,  Table S3 and Table S4. S20  Table S4. To illustrate the accuracy of the non-Faradaic capacitancecurrent model in fitting to the background current of the experimental data, Figure S9 shows an overlay of a pH 6.0 experimental dataset and the corresponding simulation of the capacitance using the best-fit parameter values.

Figure S9
(Red) The current trace for pH 6.0 experiment 3, and (blue) the corresponding simulation of the capacitance using the best-fit parameter values reported in Table S4.    Figure 4 demonstrates that clear differences in FTACV harmonics are observed in the modelling of very fast electrode kinetics, which mimic the responses predicted in the reversible Nernstian regime, and slower regimes, associated with quasi-reversible processes. For reversible electron transfer, the higher order harmonic signals are characterized by well separated lobes, with the current dropping to zero between each lobe. In the slower quasi-reversible kinetic regime, a loss of resolution is observed, which is particularly clear at the outermost lobe ( Figure 4). Table S3, the best fit parameters obtained for the concerted two-electron transfer (n3=2) model fall into the slow kinetic regime. Conversely, as visually represented in Figure S10, the experimental data falls into the fast kinetic regime, with the best fits of the n3=2 model failing to accurately describe the current traces of harmonics 6, 7, and 8. This is supportive of the conclusion that the n3=2 model is mis-specifed for the H2ase-MFHypD data at all pH's studied.

Figure S10
The 6th (top), 7th (middle), and 8th (bottom) harmonics for experiment 3 at (a) pH 4.0, (b) pH 5.0, (c) pH 5.0, (d) pH 6.0, (e) pH 8.0, and (f) pH 9.0 along with harmonics generated from the best fit parameters of the concerted two-electron transfer model to the experimental data reported in Table S3. It is common practice to verify that the faradaic current arises from protein film formation, and is not an artefact, by undertaking voltametric experiments with 'blank' electrodes before and after adhering the protein. The total (AC plus DC) current voltammograms derived from experiments conducted at pH 9.0 before and after protein film formation is shown in Figure S11; this verifies the protein had adhered and the observed faradaic signal results from the protein.

Inference and MCMC for Eapp versus pH
To probe the relationship between the separation of app1 0 and app2 0 with pH we have used the infrastructure of the 'PINTS' Python library to fit the pH vs Eapp values inferred in this work (Table 1) to Equation 2 and Equation 3 (reproduced here as Equation S34 and Equation S35 recpectively), taken from the literature (Hirst, 2006). The pKa/Ka values of these equations are defined in Figure S12, as are the true reversible potentials, 1 0 and 2 0 . When fitting the standard deviation of the noise on the reversible potentials 1 0 and 2 0 were also inferred, as a Gaussian log likelihood objective function was used. Figure S12 Square scheme from Figure 1(b) updated to define the appropriate pKa and E 0 symbols associated with each protonation and redox reaction, respectively, used to define the voltammetry of H2ase-MFHypD. This is a reproduction of Figure6(a) without arrows to indicated suggested routes across the square scheme.
As with the other inverse problems in this work, the infrastructure provide by the PINTS Python library was used and the model was written in a 'Pints.ForwardModel' class . (Clerx et al., 2019) Twenty optimizations were then performed using the C-MAES optimiser; and a uniform prior transformed by 'Pints.RectangularBoundriesTranformation' to aid fitting consistening of the following bounds: all pKa's were bound between 0.0 and 14, the true reversible potentials, 1 0 and 2 0 , were bound between -0.3 V and -0.8 V, and the standard deviation of the noise on the reversible potentials, 1 0 and 2 0 , were bound between 10 −2 and10 −8 . For an objective function the following Gaussian log likelihood function, as implemented in the 'PINTS' Python library was used: where is a vector containing the model input parameters, is a vector of the standard deviation of the noise on the model output parameters ( ), is a vector of the experimental data ( ) at each pH point, is the number of pH/data points, 0 is the number of output of the system of equations being analyzed (2 in our case), and ( ) is the simulated data at pH and input parameters . The best fit from this fitting process is given in Figure 6(b) and gave an objective function score of 149 to 3 S.F.
Posterior distributions were subsequently obtained using an adaptive Metropolis MCMC algorithm, given as algorithm 4 in (Andrieu and Thoms, 2008) and described in (Haario et al., 2001), as implemented in the 'PINTS' Python library as 'pints.HaarioACMC' (Clerx et al., 2019). Three chains were used starting at 1.15, 1.1, and 0.9 times the best fit positions previously obtained ( Figure  6(b)), chain lengths of 45000, and burn in of 20000 was determined from observation of the chain plots in Figure S13, the autocorrelation plots in Figure S14(b), and the RHAT values reported in Table S5. A uniform prior was used for all parameters covering the same ranges used for optimization and the same log-likelihood function was also used, Equation S36. The obtained posterior distributions are given in Figure S14(a). Additionally, the range of these distributions are shown in Figure 7, where they are wide enough to be shown. Pairwise distributions are shown in Figure S15.

S31
13 Insensitivity of Faradaic fitting to Uncompensated Resistance Figure S16 demonstrates the insensitivity of the faradaic current to uncompensated resistance in designated harmonics and hence our fitting method. Figure S16 (a-d, top panels) Current-potential plots of harmonics 6 to 9 of (grey line) pH 4.0 experiment 3 ("Exp") data and best fit simulated data using (red line) a concerted two-electron (n3 = 2) electron transfer model and (blue line) a concerted two-electron using best fit parameters but with an uncompensated resistance of 117Ω. (Bottom panels) Residual current-potential plots corresponding to the panel above. The parameter values used to generate the simulations are reported in Table 1 and Table S3.