Proton transport (PT), defined by the migration of a delocalized positive charge, is a ubiquitous occurrence that plays an instrumental role in numerous physiological processes in biological systems, especially in proton-coupled phosphate transporters and proton channels. Due to its crucial contribution to the functioning of these membrane proteins, a deeper understanding of how protons couple with ligand transport in the proton-coupled phosphate symporter and how proton transport differs from other ion transports - thus ensuring exceptional selectivity in the proton channel - can illuminate the protein features that may trigger enhanced or impaired transport behaviors, as well as how diseases linked to these proteins can be treated.Simulating the proton transfer (PT) process in pure solution is already a complex and nontrivial task. This problem becomes even more challenging within the context of protein systems. Molecular Dynamics (MD) simulation has emerged as a powerful tool for studying chemical reactions, especially PT in intricate protein environments. In MD simulations, PT is typically considered implicitly, as seen in simulations with empirical classical forcefields, or explicitly accounted for through hybrid Quantum Mechanics/Molecular Mechanics (QM/MM) and the Multiscale Reactive MD (MS-RMD) method. It's crucial to remember that PT can be associated with changes in a protein's conformation, ranging from brief, localized adjustments taking only several nanoseconds to more extensive, global conformational changes that can span milliseconds to minutes. The timespan of PT and its related conformational transition usually involve high free energy barriers, and as a result, can easily go beyond the limit of the vanilla unbiased MD approaches. For sufficient sampling of such rare events, enhanced sampling techniques are often required to accelerate key collective motions relevant to PT processes. The free energy profile obtained by enhanced sampling methods offers the thermodynamics property of the studied system in principle. Additional kinetic information can be inferred based on either a transition state theory assumption, a Markov state model, or a bottom-up kinetic network model. Here, we employed multiscale simulations in conjunction with enhanced sampling approaches to investigate the functional cycle of a proton-coupled phosphate symporter. We specifically focused on confirming the proposed functional role of a key residue, D324. By developing multiple possible reaction pathways within a bottom-up kinetic network model, we were able to account for different forms of phosphate. This approach allowed us to analyze the interaction between proton and phosphate gradients within each elementary reaction while making a prediction of the optimal pH condition for the phosphate uptake activity of this transporter under various phosphate concentrations. A critical part of our research involved examining the treatment of transition rate theory in the proton transfer (PT) process. We scrutinized this across multiple protein systems where one-dimensional and two-dimensional free energy profiles were available. Consequently, we derived and analyzed the source of error in each treatment. Moreover, we used MS-RMD simulations to predict the pH-dependent conductance behavior of a proton channel. The quantitative permeability analysis for proton, tetramethylammonium (TMA) cation and CH3SO3- anion proposed a perspective on the origin of the proton selectivity of the WT channel and the anion selectivity of the mutated channel D112N.