Setup of System I: The crystallographic structure of vibrio harveyi luciferase (PDB ID: 3FGC) was utilized for model construction. The structure comprises two heterodimers, 1 and 2, each containing two homologous subunits, α and β, with flavin mononucleotide (FMN) binding to the α-subunit of heterodimer 1. For the sake of computational efficiency, only heterodimer 1, including chains A and B, was retained in the model. The missing residues (residue IDs: 284-289) in chain A were added using the SWISS-MODEL online server (available at: swissmodel.expasy.org). Since the missing residues (residue IDs: 319-332) in chain B are located at the end of the chain and are far away from the fluorophore, they were reasonably ignored. The FMN in the crystal structure was modeled as HFOH. The protonation states of titratable residues (histidine, glutamic acid and aspartic acid) were determined using the online PDB2PQR program (available at: server.poissonboltzmann.org) in combination with the careful visual inspection of local hydrogen-bonded networks. All glutamic acid and aspartic acid residues were deprotonated. Histidine 45 in chain A and histidine 45 in chain B were doubly protonated, while the histidine residues 61, 82, 199, 215, 224, 234, and 249 in chain A and residues 29, 81, 82, 132, 145, 161, 209, and 224 in chain B were protonated at the ε position, histidine residues 44, 150, and 285 in chain A and residues 76 and 215 in chain B were protonated at the δ position. The Amber ff19SB force field was employed for canonical protein residues, while the general amber force field (GAFF) was used to describe HFOH. Charge parameters of HFOH were obtained from the RESP method at the HF/6-31G* level with Gaussian 16. The parmchk2 utility from Am-berTools was used to generate missing parameters. Sodium ions were added to the protein surface to neutralize the overall charge of the systems. The resulted system was solvated in a rectangular box of TIP3P waters extending up to a minimum cutoff of 12 Å from the protein surface.