The process of proton solvation and transport is involved in a wide array of phenomena. Although seemingly simple, it is a complex process that continues to engage research efforts. One unique aspect of the molecular mechanism of proton solvation and transport is that the hydrated proton is able to hop between neighboring water molecules by “Grotthuss shuttling”, giving rise to anomalously high diffusion rate of protons compared to other simple cations. In addition, the hydrated proton is able to dynamically delocalize the excess charge defect over multiple water molecules within several solvation shells. Inside proteins, the structure and dynamics of the residues and water molecules that mediate proton transport are strongly influenced by the electrostatic and hydrophobic environment in the protein cavities. Hence, proton transport in proteins is very different from that in aqueous bulk solution and displays many interesting behaviors. Explicitly simulating proton transport in proteins is an inherently challenging problem. It requires both the explicit treatment of the excess proton, including its charge defect delocalization and Grotthuss shuttling through inhomogeneous moieties (water and amino acid residues), and extensive sampling of slow conformational changes of both the protein and water clusters inside protein cavities. In recent years, the self-consistent charge density functional tight binding (SCC-DFTB) method has been increasingly applied to study proton transport (PT) in biological environments. However, recent studies revealing some significant limitations of SCC-DFTB for proton and hydroxide solvation and transport in bulk aqueous systems call into question its accuracy for simulating PT in biological systems. The current work benchmarks the SCC-DFTB/MM method against more accurate DFT/MM by simulating PT in a synthetic leucine-serine channel (LS2), which emulates the structure and function of biomolecular proton channels. It is observed that SCC-DFTB/MM produces over-coordinated and less structured pore water, an over-coordinated excess proton, weak hydrogen bonds around the excess proton charge defect and qualitatively different PT dynamics. Similar issues are demonstrated for PT in a carbon nanotube, indicating that the inaccuracies found for SCC-DFTB are not due to the point charge based QM/MM electrostatic coupling scheme, but rather to the approximations of the semiempirical method itself. The results presented in this work highlight the limitations of the present form of the SCC-DFTB/MM approach for simulating PT processes in biological protein or channel-like environments, while providing benchmark results that may lead to an improvement of the underlying method. An alternative approach that explicitly accounts for the reactive nature of the hydrated excess proton is multiscale reactive molecular dynamics (MS-RMD) method. In this approach, quantum mechanical forces from targeted quantum mechanics/molecular mechanics (QM/MM) calculations are bridged, in a multiscale fashion via a variational mathematical framework, into the reactive MD algorithm for the dynamics of system nuclei, thus allowing efficient and accurate description of Grotthuss shuttling and charge delocalization during PT. Herein, we have used a synthesis of the MS-RMD, QM/MM and classical MD simulations to study the PT mechanism in the influenza A virus M2 channel (AM2), which is crucial in the viral life cycle. Despite many previous experimental and computational studies, the mechanism of the activating process in which proton permeation acidifies the virion to release the viral RNA and core proteins is not well understood. We report, to our knowledge, the first complete free-energy profiles for PT through the entire AM2 transmembrane domain at various pH values, including explicit treatment of excess proton charge delocalization and shuttling through the His37 tetrad. The free-energy profiles reveal that the excess proton must overcome a large free- energy barrier to diffuse to the His37 tetrad, where it is stabilized in a deep minimum reflecting the delocalization of the excess charge among the histidines and the cost of shuttling the proton past them. At lower pH values the His37 tetrad has a larger total charge that increases the channel width, hydration, and solvent dynamics, in agreement with recent 2D-IR spectroscopic studies. The PT barrier becomes smaller, despite the increased charge repulsion, due to backbone expansion and the more dynamic pore water molecules. The calculated conductances are in quantitative agreement with recent experimental measurements. In addition, the free-energy profiles and conductances for PT in several mutants provide insights for explaining our findings and those of previous experimental mutagenesis studies. Another important challenge for simulation of PT in biomolecular systems is a quantitative description of the protonation and deprotonation process of amino acid residues. Despite the seeming simplicity of adding or removing a positively charged hydrogen nucleus, simulating the actual protonation/deprotonation process is inherently difficult. In a recent paper (J. Chem. Theory Comput. 2014, 10, 2729−2737), a multiscale approach was developed to map high-level quantum mechanics/molecular mechanics (QM/MM) data into a MS-RMD model in order to describe amino acid deprotonation in bulk water. Herein, the fitting approach (called FitRMD) was extended to create MS-RMD models for ionizable amino acids within proteins. The resulting models are shown to faithfully reproduce the free energy profiles of the reference QM/MM Hamiltonian for PT inside an example protein, the ClC-ec1 H+/Cl− antiporter. Moreover, we show that the resulting MS-RMD models are computationally efficient enough to then characterize more complex 2-dimensional free energy surfaces due to slow degrees of freedom such as water hydration of internal protein cavities that can be inherently coupled to the excess proton charge translocation. The FitRMD method is thus shown to be an effective way to map ab initio level accuracy into a much more computationally efficient reactive MD method in order to explicitly simulate and quantitatively describe amino acid protonation/deprotonation in proteins. The combination of FitRMD and our other multiscale methods is herein used to study proton transport in aa3-type cytochrome c oxidase (CcO), which is the terminal enzyme in the respiratory electron transfer chain in the inner membrane of mitochondria and plasma membrane of bacteria. CcO catalyzes the reduction of O2 to H2O and couples the free energy of this exergonic reaction to the pumping of protons across the membrane, creating a transmembrane proton electrochemical gradient that drives, for example, ATP synthesis. We have used MS-RMD simulations to explicitly characterize (with free energy profiles and calculated rates) the internal PT events that enable pumping and chemistry during the A→PR→F transition in the aa3-type CcO. Our results show that PT from amino acid residue E286 to both the pump loading site (PLS) and to the binuclear center (BNC) are thermodynamically driven by electron transfer from heme a to the BNC, but that the former (i.e., pumping) is kinetically favored while the latter (i.e., transfer of the chemical proton) is rate-limiting. The calculated rates are in quantitative agreement with experimental measurements. The back flow of the pumped proton from the PLS to E286 and from E286 to the inner side of membrane are prevented by the fast reprotonation of E286 through the D-channel and large free energy barriers for the back flow reactions. PT from E286 to the PLS through the hydrophobic cavity (HC) and from D132 to E286 through the D-channel are found to be strongly coupled to dynamical hydration changes in the corresponding pathways. This work presents a comprehensive description of the key steps in the proton pumping mechanism in CcO. Lastly, we have used a synthesis of the MS-RMD, QM/MM and classical MD simulations to explicitly characterize (with free energy profiles and calculated conductances) the activation mechanism of influenza M2 channel from high to intermediate pH condition. Starting from a recently resolved crystal structure in high pH condition (4QKL), we found that the Trp41 side chains in the high pH condition block the channel, and the channel is gradually activated as the pH in the viral exterior is lowered, as a result of the increased positive charges on the His37 tetrad. The qualitative asymmetry of the PT free energy profile in high pH condition explains why the inward proton flux is allowed when the pH is low in viral exterior and high in viral interior, but outward proton flux is prohibited when the pH gradient is reversed. Our results also demonstrate that the amphipathic helices do not change the proton conduction mechanism in the AM2 transmembrane domain significantly.