Protein residue ionization often plays a vital role in protein function, such as enzyme catalysis and pH-dependent structural changes. Inadequate conformational and protonation state sampling in molecular dynamics (MD) simulations has limited the accuracy of predictions of residue pKa values and structural change concomitant to residue ionization state change. As part of an ongoing effort in the Yang lab to address this conformational sampling issue, we developed a novel enhanced sampling MD simulation method, called adaptive dynamic reporting (ADR). In this thesis, we show that the ADR method has excellent accuracy, precision, and convergence time for a wide range of simulation parameter sets and order parameters, compared to previously established methods. Next, using ADR-based high-order orthogonal space tempering (OST) methods, we studied a number of ionization-dependent properties for all human oxidized C62A/C69A/C73A thioredoxin ionizable residues. To carry our this study, we developed a novel alchemical OST pKa titration scheme: 1) to sample the relevant protonation microstates of thioredoxin, and 2) to calculate the microscopic pKa values, protonation microstate populations, titration curves, and apparent pKa for all human oxidized C62A/C69A/C73A thioredoxin ionizable residues. The OST pKa simulations accurately predict the apparent pKa values for each of these residues. Moreover, from our OST pKa simulations we identified novel thioredoxin conformations and dynamics corresponding to biologically essential ionization states. These long-timescale, ionization-coupled motions observed were in good agreement with NMR spectroscopy observations of homologous thioredoxin. In addition, we were able to utilize the observed the ionization-coupled dynamics of human oxidized C62A/C69A/C73A thioredoxin to identify possible sources of the pKa shifts for several titratable residues. We then studied ionization coupling between residues using their deprotonation free energies. Several interacting residues showed nonlinear deprotonation free energies over various pH ranges. In summary, there are four main contributions in this dissertation. First, we developed and analyzed a novel, rigorous, and widely-applicable enhanced sampling method, ADR. Second, we developed a novel titration scheme, and in combination with OST pKa simulations, calculated and analyzed all of the residue microscopic pKa values in human oxidized C62A/C69A/C73A thioredoxin. Third, we established how the hierarchical dynamics of thioredoxin are constructed to respond to external pH perturbations by characterizing novel protonation-dependent dynamics of thioredoxin. Finally, we extensively studied inter-residue ionization cooperativity in thioredoxin by calculating residue deprotonation free energies and by analyzing the dynamics of thioredoxin coupled to the ionization of these residues.