**3. Discussion**

The absence of the relevant UCP structure hinders the understanding of UCP2 transport mechanisms and biological functions. To find a possible remedy for the current situation, we turned to homology modeling using ANT crystallographic structure and experimental UCP2 NMR structure as starting structures for long microsecond MD simulations of UCP2 proteins (UCP2h and UCP2NMR). We showed by microsecond MD simulations that UCP2h structure is almost water impermeable (*Pf* = 2.0 × 10−<sup>16</sup> cm<sup>3</sup> s<sup>−</sup>1), with a water osmotic permeability coefficient lower by three orders of magnitude than the corresponding UCP2NMR structure after 2 μs of MD simulation time (*Pf* = 1.3 × 10−<sup>13</sup> cm<sup>3</sup> s<sup>−</sup>1). This is a consequence of the salt bridge network formation at the matrix side of the protein in UCP2h, formed by arginine/lysine and aspartate/glutamate residues, as well as conservation of the threefold symmetry between negatively charged residues at the matrix side (EG motif),

similar to water-impermeable ANT proteins [44,56,57]. In contrast, the stable formation of the identical salt bridges and the EG motif at the matrix side of UCP2NMR structure were not observed, thus enabling water to almost freely diffuse across the protein. Finally, we showed that ATP binds to all three arginine residues in the UCP2 cavity only in the case of UCP2h structure, in agreemen<sup>t</sup> with previous experimental observations [53,54]. We should mention here that MD simulations of UCP2 proteins are not fully converged (which is almost impossible with the current computational power) and that UCP2h structure is possibly not the global minimum. Still, we believe that the presented UCP2h structure is relevant for future studies since it corresponds well to available experimental data and that its behavior is described sufficiently well.

The described UCP2h structure can be easily transferred to membranes formed by DOPC, DOPE and cardiolipin, which would then more realistically describe the lipid composition of inner mitochondrial membranes and our model membranes. However, at this stage, we believe that further complication of a model system is not necessary, as we have shown previously in MD simulations of ANT proteins in different, more complex environments [44]. Due to a high complexity of membrane protein/lipid bilayer systems and a necessity for exceedingly long simulation times, pitfalls resulting from incomplete coverage of whole phase space in heterogeneous systems, even at the microsecond timescales and relatively small bilayers, are most likely.

To test the validity of MD simulation results regarding UCP2 functionality, we compared the binding of UCP2 activator AA to both models, UCP2h and UCP2NMR. Additionally, we measured the activation of the recombinant UCP2 due to the AA binding. A larger membrane conductance of the UCP2 in the presence of AA compared to the UCP2-R60S, in which the putative AA− binding site was modified, is in full agreemen<sup>t</sup> with MD simulation results. It convincingly shows that the suggested AA− binding region around R60 is available only in the UCP2h protein structure and that in UCP2-R60S the binding of AA− is prevented. This shows that the UCP2h structure correctly predicts the binding of AA− and is sensitive to selective mutations, such as R60S. In contrast, the similar binding of FA to the UCP2NMR structure, both the experimental one and after 2 μs of MD simulation time is not optimal, thus further questioning the stability and physiological relevance of experimental NMR structure.

The proposed MD model of UCP2h can largely contribute to the investigation of the UCP2-mediated proton transport mechanism in general and the protein binding site for FA in particular. Currently, two basic mechanisms are proposed for UCP1-mediated proton transport (see Introduction). A crucial difference between the "FA cycling" and "FA shuttle" hypotheses is a localization of the binding site for the activating FA. The "FA shuttle" model states that FA− binds inside the pore from the cytosolic side of UCP1 and transfers protons by shuttling from the cytosolic to the matrix side of mitochondria [58]. In contrast, the "FA cycling" hypothesis suggests that FA− binds to protein from the matrix side and is transported to the intermembrane space in its deprotonated form. Our initial results that show that FA can bind from the matrix site contradict the results of patchclamp experiments on mitoplasts, showing the inability of UCP1 to bind FAs on the matrix side [18]. Further developments of the UCP2 model will pave the way towards clarifying this question and the complete transport mechanism of UCP2.

In conclusion, we propose the protocol involving homology modelling of UCP2 protein based on ANT structure and subsequent microsecond molecular dynamics simulations. It yields a functionally relevant structure, which can be used for future mechanistic studies of proton/FA transfer in mitochondria.

#### **4. Materials and Methods**

## *4.1. Simulation Details*

To analyze the UCP2 protein structure, schematically represented in Figure S1, we simulated and compared two different structures of UCP2 protein and bovine ANT protein as a reference. The simulated structures include (a) the published NMR structure of mUCP2 (PDB code 2LCK, organism *Mus musculus*) [19] determined by NMR molecular fragment replacement, (b) the homology modeled structure of the primary sequence of UCP2 protein to crystallographic structure of ANT protein (PDB code 1OKC, 2.2 Å, high resolution, organism *Bos taurus*) [38] whose homology to mUCP2 is 24%, and (c) a crystallographic structure of ANT protein, which served as a template for the homology modeling of UCP2. Moreover, we used it as a reference structure for investigated UCP2 structures.

Three different systems were prepared using the membrane builder module of CHARMM-GUI (http://www.charmm-gui.org/) [59–61] from three UCP2/ANT starting protein structures immersed in 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) membrane containing 230 lipid molecules, 28,750 water molecules, and 15/19 chloride ions to neutralize the net charge of UCP2/ANT protein. DOPC was selected as it represents one of the main lipids of inner mitochondrial membranes. The homology modeled structure of UCP2 and missing residues (compared to the crystallographically determined structures of ANT)—residue 1 and residues 294–297 for 1OKC and residues 1–13 for 2LCK—were added using USCF Chimera program [62]. All arginines and lysines were prepared in their protonated forms, histidines and cysteines in their neutral forms, with glutamates and aspartates being in their deprotonated forms. CHARMM-GUI membrane builder minimization and equilibration procedure was used for all systems [59]. After equilibration, we simulated each system for 2.0 μs using unbiased all-atom molecular dynamics (MD) simulations in a periodic rectangular box of 9.5 nm × 9.5 nm × 13.5 nm with a time step of 2 fs with CHARMM36m force field [63] and TIP3P water model [64].

In order to check the potential binding site of the fatty acid anion in UCP2, we also simulated an additional four systems of UCP2 protein structures with added arachidonic acid anion (AA−). AA was used since it shows a substantial effect on the proton conductance across mitochondrial membranes [10]. AA− was added to the binding site arginine R60 in three systems: (a) the homology modelled structure of UCP2 (UCP2h) simulated after 2 μs, (b) the NMR structure after initial equilibration which closely resembles experimental NMR structure, and (c) the NMR experimental structure (UCP2NMR) simulated after 2 μs. In the fourth system, AA− was added to the mutated binding site (R60 to S60) in UCP2h simulated after 2 μs (UCP2-R60h). Each of the four AA− containing systems mentioned above was simulated for 500 ns and analyzed further.

All production MD simulations were performed in the isobaric-isothermal ensemble (NPT) at *T* = 310 K, which was maintained via Nosé−Hoover thermostat [65,66] independently for three groups: DOPC, water/ions, and protein subsystem with a coupling constant of 1.0 ps<sup>−</sup>1. The pressure was set to 1.013 bar and controlled with semi-isotropic Parrinello −Rahman barostat [67] with a time constant for pressure coupling of 5 ps<sup>−</sup>1. Periodic boundary conditions (PBC) were imposed in all three directions, with long range electrostatic interactions calculated by the particle-mesh Ewald (PME) method [68] with real space Coulomb interactions cut-off at 1.2 nm using a Fourier spacing of 0.12 nm and Verlet cut-off scheme. All simulations were propagated with GROMACS 2018.4 software package [69] and visualized with VMD (Visualize Molecular Dynamics) program [70].
