Overview

---

The metabolic processes in living organisms inherently possess stochastic elements and intrinsic noise, rendering the attainment of a strictly analytical or precisely deterministic solution a formidable challenge in the realm of biomathematics. In response to this complexity, biologists and mathematicians frequently employ kinetics and statistics as powerful tools to establish mathematical models that encapsulate the dynamics of biochemical reactions ^[1]. In the context of elucidating the operational intricacies of the formaldehyde metabolism pathways and NCD synthesis within E. coli, and, more significantly, in comprehending the ramifications of alterations in the efficiency and activity of the enzymes governing these processes on the system's overall behaviour, we advocate for an in-depth, bottom-up approach. A kinetic modelling framework, grounded in the metabolic pathway, is essential for scrutinizing this dynamic system.

Biochemical Description and Analysis

---

In the examination of the NCD-dependent formaldehyde metabolic system, it is imperative to adopt a biochemist's perspective. Complex biochemical networks are conventionally employed to represent intricate metabolic pathways, and herein, we provide a structural formula-based representation of the NCD-dependent formaldehyde pathway (Fig. 1). This figurative framework enables us to dissect the system's progression systematically. The influx of formaldehyde into the cellular matrix is driven by concentration gradients, and, guided by our reengineered metabolic pathway, it undergoes a transformation within the cell, effectively functioning as a formaldehyde processing conduit, akin to a specialized pipeline.

Biochemical description

---

In a biochemical context, the progression of formaldehyde metabolism unfolds as follows:

1. Exogenous Formaldehyde Passive Diffusion: The initial stage entails the passive diffusion of formaldehyde (HCHO), a small hydrophilic molecule, across the cell membrane driven by concentration gradients. This diffusion process is reversible and is distinguished by subscripts ex and in denoting the formaldehyde's presence in distinct cellular compartments.

   $\mathrm{Formaldehyde}{\phantom{A}}_{\mathrm{ex}}\stackrel{-\rightharpoonup }{\leftharpoondown -}\mathrm{Formaldehyde}{\phantom{A}}_{\mathrm{in}}$$\ce{Formaldehyde_{ex} <=> Formaldehyde_{in}}$

2. Oxidation of Cellular Formaldehyde: Following diffusion, formaldehyde encounters oxidation catalyzed by formaldehyde dehydrogenase (FADH). This catalytic reaction results in the conversion of hazardous formaldehyde into formate. Notably, the non-natural cofactor NCD, an artificial analog of NAD, plays a pivotal role as an electroreceptor in facilitating this reaction. FADH catalyzes this reaction, which can be represented as a reversible process.

   $\mathrm{Formaldehyde}+\mathrm{NCD}{\phantom{A}}^{+}+\mathrm{H}{\phantom{A}}_2\mathrm{O}\stackrel{\mathrm{FADH}}{\stackrel{-\rightharpoonup }{\leftharpoondown -}}\mathrm{Formate}+\mathrm{NCDH}+\mathrm{H}{\phantom{A}}^{+}$$\ce{Formaldehyde + NCD+ + H2O <=>[FADH] Formate + NCDH + H+}$

   A structural comparison between NAD (Nicotinamide Adenine Dinucleotide) and NCD (Nicotinamide Cytosine Dinucleotide) is depicted in Fig. 2, revealing the distinctions marked by cyan and tawny rectangles. While NCD and NAD share a common dinucleotide backbone covalently linked by pyrophosphate bonds, their different functional nitrogenous bases impart varying biochemical roles in vivo ^[2].

3. Oxidation of Formate: Formate, the product of formaldehyde oxidation, undergoes further transformation into carbon dioxide, with this reaction also being dependent on NCD. This process involves an artificially mutated formate dehydrogenase (FDH*), distinguished by an asterisk, indicating differences from the wild-type FDH.

   $\mathrm{Formate}+\mathrm{NCD}{\phantom{A}}^{+}\stackrel{\mathrm{FDH}{\phantom{A}}^{\ast }}{\stackrel{-\rightharpoonup }{\leftharpoondown -}}\mathrm{CO}{\phantom{A}}_2+\mathrm{NCDH}+\mathrm{H}{\phantom{A}}^{+}$$\ce{Formate + NCD+ <=>[FDH^{*}] CO2 + NCDH + H+}$

4. Incorporating Downstream Products: The challenge of incorporating the downstream product of formaldehyde into the existing natural metabolic pathways is addressed by guiding carbon dioxide (CO2) towards the Tricarboxylic Acid Cycle (TCA) through a series of steps facilitated by a mutant-type malic enzyme (ME*). By bibliographic retrieval, we find that this coupled reaction can be divided into two oxaloacetate-involved reactions ^[2,3,4]. This transformation of CO2 into malate is achieved through a progression known as the "Hydrocarboxylation of Pyruvate."

   $\mathrm{Pyruvate}+\mathrm{CO}{\phantom{A}}_2+\mathrm{NCDH}+\mathrm{H}{\phantom{A}}^{+}\stackrel{\mathrm{ME}{\phantom{A}}^{\ast }}{\to }\mathrm{Malate}+\mathrm{NCD}{\phantom{A}}^{+}$$\ce{Pyruvate + CO2 + NCDH + H+ ->[ME^*] Malate + NCD+}$

   This complex reaction comprises two distinct sub-reactions:

   a. Carboxylation of Pyruvate: In the first sub-reaction, pyruvate is carboxylated, yielding oxaloacetate intermediates. This step enables the immobilization of carbon dioxide.

   $\mathrm{Pyruvate}+\mathrm{CO}{\phantom{A}}_2\stackrel{\mathrm{ME}{\phantom{A}}^{\ast }}{\to }\mathrm{Oxaloacetate}$$\ce{Pyruvate + CO2 ->[ME^*] Oxaloacetate}$

   b. Reduction of Oxaloacetate Intermediate: In the second sub-reaction, the oxaloacetate intermediate is reduced, with the artificial proton donor NCDH facilitating the reaction and regenerating the oxidized form of NCD. The primary product, malate, participates in various biochemical processes in vivo, including the TCA cycle. Significantly, this redesigned formaldehyde metabolic pathway presents an eco-friendly approach for malate production, reducing the reliance on hydrogen (H₂) as seen in conventional industrial processes ^[4].

   $\mathrm{Oxaloacetate}+\mathrm{NCDH}+\mathrm{H}{\phantom{A}}^{+}\stackrel{\mathrm{ME}{\phantom{A}}^{\ast }}{\to }\mathrm{Malate}+\mathrm{NCD}{\phantom{A}}^{+}$$\ce{Oxaloacetate + NCDH + H+ ->[ME^*] Malate + NCD+}$

The mechanism of hydrocarboxlyation reaction

---

In the elucidation of the hydrocarboxylation reaction mechanism, we present a rational model based on a combination of structural analysis and comprehensive bibliographic research ^[2,3,4,5]. The process commences with pyruvate expelling a proton under the influence of a Brønsted base, leading to the conjugation of the lone pair of the γ-carbon with the β-carbonyl group. Employing resonance structures, the enolate anion isomer arises from the carbonyl anion, offering enhanced suitability for nucleophilic attacks. Subsequently, the enolate form of the pyruvate anion participates in a nucleophilic addition with carbon dioxide originating from the formaldehyde oxidation pathway, resulting in the formation of the oxaloacetate anion. This oxaloacetate anion then undergoes further reduction to yield the malate anion, all facilitated by the enzymatic action of ME^*. Finally, the malate anion is protonated to attain a stable species, culminating in the hydrocarboxylation reaction (Fig. 3).

Chemical Kinetic Model

---

Model establishment and analysis

---

Kinetic analysis stands as a time-honoured and highly effective approach for the anticipation and computational simulation of intricate biochemical networks. In the context of our study, the rate variables within this system, e.g., formaldehyde, NCD, CO₂, are no longer represented directly but have instead been replaced by a combination of concentrations for specific chemical species. This transformation allows us to embark on the crucial task of simulating the dynamic behaviour of the newly redesigned metabolic network (Fig. 4).

For a general reaction, $s$$s$ substrates with concentration $[S{]}_1,[S{]}_2,[S{]}_i,...,[S{]}_s$$[S]_1,[S]_2,[S]_i,...,[S]_s$ are converted into $p$$p$ products with concentration $[P{]}_1,[P{]}_2,[P{]}_j,...,[P{]}_p$$[P]_1,[P]_2,[P]_j,...,[P]_p$, we have

$\mathrm{S}{\phantom{A}}_1+\mathrm{S}{\phantom{A}}_2+\dots +\mathrm{S}{\phantom{A}}_{\mathit{s}}\underset{{\mathit{v}}_{\mathit{-}}}{\stackrel{{\mathit{v}}_{\mathit{+}}}{\stackrel{-\rightharpoonup }{\leftharpoondown -}}}\mathrm{P}{\phantom{A}}_1+\mathrm{P}{\phantom{A}}_2+\dots +\mathrm{P}{\phantom{A}}_{\mathit{p}}$$\ce{S_1 + S_2 + ... + S_{\mathit{s}} <=>[\mathit{v_{+}}][\mathit{v_{-}}] P_1 + P_2 + ... + P_{\mathit{p}}}$

$v=v_{+}-v_{-}=k_{+}\prod _{i=1}^s{S}_i^{n_i}-k_{-}\prod _{j=1}^p{P}_j^{n_j},$$v=v_{+}-v_{-}=k_{+}\prod_{i=1}^{s}S^{n_i}_i - k_{-}\prod_{j=1}^{p}P^{n_j}_j,$

where $n_i$$n_i$ and $n_j$$n_j$ define molecularities of substrate $S_i$$S_i$ and product $P_j$$P_j$ respectively. We can describe the dynamics of the HCHO metabolic pathway by the linear combination of rates, here is a set of ordinary differential equations assigned to particular species as follows.

$\frac{\mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{ex}}}{\mathrm{d}t}=-v_{+}^{(0)}+v_{-}^{(0)}$$\frac{\mathrm{d\left[Formaldehyde\right]_{ex}}}{\mathrm{d}t}=-v^{(0)}_{+}+v^{(0)}_{-}$

$\frac{\mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{in}}}{\mathrm{d}t}=v_{+}^{(0)}-v_{-}^{(0)}-v_{+}^{(1)}+v_{-}^{(1)}$$\frac{\mathrm{d\left[Formaldehyde\right]_{in}}}{\mathrm{d}t}=v^{(0)}_{+}-v^{(0)}_{-}-v^{(1)}_{+}+v^{(1)}_{-}$

$\frac{\mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{in}}}{\mathrm{d}t}=v_{+}^{(0)}-v_{-}^{(0)}-v_{+}^{(1)}+v_{-}^{(1)}$$\frac{\mathrm{d\left[Formaldehyde\right]_{in}}}{\mathrm{d}t}=v^{(0)}_{+}-v^{(0)}_{-}-v^{(1)}_{+}+v^{(1)}_{-}$

$\frac{\mathrm{d}\left[{\mathrm{CO}}_2\right]}{\mathrm{d}t}=v_{+}^{(2)}-v_{-}^{(2)}-v^{(3)}$$\frac{\mathrm{d\left[CO_{2}\right]}}{\mathrm{d}t}=v^{(2)}_{+}-v^{(2)}_{-}-v^{(3)}$

$\frac{\mathrm{d}\left[{\mathrm{CO}}_2\right]}{\mathrm{d}t}=v_{+}^{(2)}-v_{-}^{(2)}-v^{(3)}$$\frac{\mathrm{d\left[CO_{2}\right]}}{\mathrm{d}t}=v^{(2)}_{+}-v^{(2)}_{-}-v^{(3)}$

$\frac{\mathrm{d}\left[{\mathrm{NCD}}^{+}\right]}{\mathrm{d}t}=-\frac{\mathrm{d}\left[\mathrm{NCDH}\right]}{\mathrm{d}t}=-v_{+}^{(1)}+v_{-}^{(1)}-v_{+}^{(2)}+v_{-}^{(2)}+v^{(3)}$$\frac{\mathrm{d\left[NCD^{+}\right]}}{\mathrm{d}t}=-\frac{\mathrm{d\left[NCDH\right]}}{\mathrm{d}t}=-v^{(1)}_{+}+v^{(1)}_{-}-v^{(2)}_{+}+v^{(2)}_{-}+v^{(3)}$

$\frac{\mathrm{d}\left[\mathrm{Malate}\right]}{\mathrm{d}t}=v^{(3)}-v^{(4)}+v^{(5)}$$\frac{\mathrm{d\left[Malate\right]}}{\mathrm{d}t}=v^{(3)}-v^{(4)}+v^{(5)}$

Assign law of mass action to the above system of ODE, we have

$v_{+}^{(0)}=k_{+}^{(0)}[{\mathrm{Formaldehyde}}_{\mathrm{ex}}]$$v^{(0)}_{+}=k^{(0)}_{+}[\mathrm{Formaldehyde_{ex}}]$

$v_{-}^{(0)}=k_{-}^{(0)}[{\mathrm{Formaldehyde}}_{\mathrm{in}}]$$v^{(0)}_{-}=k^{(0)}_{-}[\mathrm{Formaldehyde_{in}}]$

$v_{+}^{(1)}=k_1[{\mathrm{Formaldehyde}}_{\mathrm{in}}][\mathrm{NCD}]$$v^{(1)}_{+}=k_1[\mathrm{Formaldehyde_{in}}][\mathrm{NCD}]$

$v_{-}^{(1)}=k_{-1}[\mathrm{Formate}][\mathrm{NCDH}]$$v^{(1)}_{-}=k_{-1}[\mathrm{Formate}][\mathrm{NCDH}]$

$v_{+}^{(2)}=k_2[\mathrm{Formate}][\mathrm{NCD}]$$v^{(2)}_{+}=k_2[\mathrm{Formate}][\mathrm{NCD}]$

$v_{-}^{(2)}=k_{-2}[{\mathrm{CO}}_2][\mathrm{NCDH}]$$v^{(2)}_{-}=k_{-2}[\mathrm{CO_2}][\mathrm{NCDH}]$

$v^{(3)}=k_3[{\mathrm{CO}}_2][\mathrm{Pyruvate}][\mathrm{NCDH}],$$v^{(3)}=k_3[\mathrm{CO_2}][\mathrm{Pyruvate}][\mathrm{NCDH}],$

where $k_{+}^{(i)}$$k^{(i)}_{+}$ and $k_{-}^{(i)}$$k^{(i)}_{-}$ denote the reaction rate constant of the forward reaction and reverse reaction respectively. To accomplish this, we capitalize on the capabilities of the SimBiology tool from the MathWorks toolbox. Notably, SimBiology is a tool commonly employed in the dynamic modelling of pharmacokinetics, an area that exhibits significant parallels with the characteristics of metabolic networks. Consequently, it becomes evident that SimBiology can yield precise and numerical outcomes in our analytical pursuits, as is further substantiated in our subsequent model optimization endeavours. This choice of methodology not only affords a comprehensive understanding of the system but also facilitates model refinement, contributing to the robustness and accuracy of our analysis (Fig. 5).

Biochemical network simulation using SimBiology

---

Fig. 6 consists of five subgraphs (a)-(e) that collectively elucidate the dynamics of various chemical species over time or other species. This metabolic system behaves in the temporal evolution of concentrations for multiple interacting species (Fig. 6a), and a notable trend emerges, where the concentration of environmental/cytoplasmic formaldehyde decreases progressively with time. This decline is juxtaposed with the behaviour of other species, which exhibit a distinctive peak-shaped pattern, except for the constitutive expressed NCD. These species show an initial increase in concentration, followed by a subsequent decrease. This unique trend underscores the intricate and dynamic nature of the chemical interactions within the system, revealing a delicate balance of processes that drive these changes over time.

In contrast, Fig. 6b-e delves into the phase diagrams, illustrating the interdependencies and fluctuations of specific species concerning the concentrations of environmental formaldehyde, cytoplasmic formaldehyde, malate, and NCD. The phase graphs offer a comprehensive overview of the relationships and trends among these chemical components, facilitating a deeper understanding of their interactions and potentially uncovering critical insights into the system under study.

Local sensitivity analysis of the biochemical network

---

The time-dependent sensitivities of environmental formaldehyde $S(\mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{ex}},k)$$S\left(\mathrm{d\left[Formaldehyde\right]_{ex}}, k\right)$ with respect to each parameter (e.g. $k_{+}^{(0)},k_{-}^{(0)}$$k^{(0)}_{+},k^{(0)}_{-}$) value are the time-dependent derivatives.

$\begin{array}{c}S(\mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{ex}},k_{+}^{(0)})=\frac{\partial \mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{ex}}}{\partial k_{+}^{(0)}}\\ S(\mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{ex}},k_{-}^{(0)})=\frac{\partial \mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{ex}}}{\partial k_{-}^{(0)}}\\ ...\\ S(\mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{ex}},k_4)=\frac{\partial \mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{ex}}}{\partial k_4},\end{array}$$S\left(\mathrm{d\left[Formaldehyde\right]_{ex}},k^{(0)}_{+}\right)=\frac{\partial\mathrm{d\left[Formaldehyde\right]_{ex}}}{\partial k^{(0)}_{+}}\\ S\left(\mathrm{d\left[Formaldehyde\right]_{ex}},k^{(0)}_{-}\right)=\frac{\partial\mathrm{d\left[Formaldehyde\right]_{ex}}}{\partial k^{(0)}_{-}}\\ ...\\ S\left(\mathrm{d\left[Formaldehyde\right]_{ex}},k_{4}\right)=\frac{\partial\mathrm{d\left[Formaldehyde\right]_{ex}}}{\partial k_{4}},$

where the numerator $\partial \mathrm{d}{\left[\mathrm{Formaldehyde}\right]}_{\mathrm{ex}}$$\partial \mathrm{d\left[Formaldehyde\right]_{ex}}$ is the output and the denominators $\partial k$$\partial k$ are the inputs to sensitivity analysis. local sensitivity analysis (LSA) was conducted using SimBiology, a platform that combines ordinary differential equation (ODE) solvers with complex-step approximation techniques. This approach allowed us to explore the system's sensitivity to parameter variations.

In conjunction with our previous analysis of dynamic trends of the metabolism system, we performed a local sensitivity analysis using SimBiology. The results revealed crucial insights into the system's behaviour. Specifically, we observed that the rate constants kf_0 and kr_0 exert a significant influence on the concentration of environmental formaldehyde, underscoring the system's dependency on diffusion progress rates (Fig. 7a). Furthermore, cellular formaldehyde concentrations exhibited high sensitivity to parameters kf_0, kr_0, and kf_ncd, highlighting the pivotal role of in vivo expression of the non-natural cofactor NCD in regulating cellular formaldehyde levels, particularly by facilitating downstream oxidation (Fig. 7c). Notably, the versatile cofactor NCD displayed heightened sensitivity to parameters kf_ncd, kf_0, kr_0 (partially), and kf_1, suggesting optimization prospects for endogenous NCD expression and the potential introduction of formaldehyde membrane transporters to enhance NCD's efficiency (Fig. 7d).

Fig. 7b shows the parameters involved in the LSA and the schematics of the external formaldehyde metabolic pathway. These LSA findings provide valuable insights for understanding and optimizing the system's behavior, offering opportunities for further research in various application domains.

In Fig. 8, we present the temporal evolution profiles of external formaldehyde, internal formaldehyde, and NCD under varying parameter conditions. To comprehensively understand the system's dynamics, we executed parameter iterations for kf_0, kr_0, kf_1, and kf_ncd. Remarkably, these iterations, performed under distinct initial conditions characterized by varying kinetic constants, led to discernable distinctions in the concentration dynamics of the species of interest. This observed variability in behaviour aligns seamlessly with the outcomes of our prior investigation into local sensitivity, as described earlier. These findings underscore the significant influence of parameter variations on the system's behaviour, offering valuable insights into its intricate dynamics and the potential for optimization.

Molecular Model

---

The aim of our project is to metabolize formaldehyde through the use of modified enzymes capable of utilizing the non-natural coenzyme NCD. Consequently, conducting molecular docking simulations with these pertinent enzymes is of paramount importance for our forthcoming wet-lab experiments. Furthermore, successful molecular docking serves as an additional validation of the effectiveness of our bioengineering efforts.

Method & Results

---

Receptor	Ligand
Ncds-2	CTP
Ncds-2	NMN

Weighted RMSD

---

Consider a molecule characterized by $N$ $N$ atoms located at positions $A={\left\{a_i\right\}}_N$ $A=\left\{a_{i}\right\}_{N}$ , with coordinates $w={\left\{w_i\right\}}_N$ $w=\left\{w_{i}\right\}_{N}$ . When presented with two groups of $N$ $N$ points, $A$ $A$ and $A^{\prime }$ $A^{\prime}$ , with respective coordinates ${\left\{a_i\right\}}_N$ $\left\{a_{i}\right\}_{N}$ and ${\left\{a_i^{\prime }\right\}}_N$ $\left\{a_{i}^{\prime}\right\}_{N}$ , which represent two different conformations of a molecule, we can establish the weighted RMSD between them ^[6] .

$\mathrm{RMSD}{(A,A^{\prime })}^2=\frac{1}{W}\sum _i{w}_i{|{\mathbf{a}}_i-{\mathbf{a}}_i^{\prime }|}^2$ $\operatorname{RMSD}\left(A, A^{\prime}\right)^{2}=\frac{1}{W} \sum_{i} w_{i}\left|\mathbf{a}_{i}-\mathbf{a}_{i}^{\prime}\right|^{2}$

The term $W$ $W$ is defined as the sum of the weights, i.e., $W=\sum _i{w}_i$ $W=\sum_{i} w_{i}$ . The statistical weights ${\left\{w_i\right\}}_N$ $\left\{w_{i}\right\}_{N}$ can be used to highlight the significance of a specific part of the structure. For instance, in the case of a protein, this could be the backbone or the side chains. More often, these weights are equivalent to the atomic masses (in this scenario, $W$ $W$ is equal to the total mass of the molecule) or they may be set to 1 (in which case, $W=N$ $W=N$ ).

We use Autodock4 for molecular docking to simulate the affinity of the enzyme for various small molecules.

Quaternion arithmetic

---

A quaternion, denoted as $Q$$Q$, can be expressed as a combination of a scalar component, $\mathrm{s}$$\mathrm{s}$, and a three-dimensional vector $\mathbf{q}={\{q_x,q_y,q_z\}}^T$$\mathbf{q}=\left\{q_{x}, q_{y}, q_{z}\right\}^{T}$, Thus, $Q$$Q$ can be written as $[s,q]$$[s, q]$. This representation is particularly useful for describing spatial transformations, especially rotations. For instance, a rotation operation can be represented by a rotation quaternion, denoted as $\hat{Q}$$\hat{Q}$, which has a unit norm. Quaternion algebra includes various operations such as multiplication, division, inversion, and norm calculation. A brief overview of quaternion arithmetic can be found in our previous work ^[7].

The rigid-body motion case

---

Here, we recapitulate the main outcome of our prior work ^[7]. We presume that a molecule A with coordinates $a_i={\{x_i,y_i,z_i\}}^T$$a_{i}=\left\{x_{i}, y_{i}, z_{i}\right\}^{T}$ is translated and rotated to new positions $A^{\prime }={\left\{{\mathbf{a}}_i^{\prime }\right\}}_N$$A^{\prime}=\left\{\mathbf{a}_{i}^{\prime}\right\}_{N}$, which are given as ${\mathbf{a}}_i^{\prime }={\mathbf{R}\mathbf{a}}_i+\mathbf{T}$$\mathbf{a}_{i}^{\prime}=\mathbf{R a}_{i}+\mathbf{T}$. Here, $\mathbf{R}$$\mathbf{R}$ is the $3\times 3$$3 \times 3$ rotation matrix and $\mathbf{T}$$\mathbf{T}$ is the translation 3 -vector. It is convenient to use the quaternion representation $Q=[s,q]$$Q=[s, q]$ of the rotation matrix $\mathbf{R}$$\mathbf{R}$. Then, the weighted RMSD between $A$$A$, the original positions, and $A^{\prime }$$A^{\prime}$, the transformed positions, can be written according to equation (Available in a centroid reference frame of C=0, Less arithmetic expressions) from Rapid determination of RMSDs corresponding to macromolecular rigid body motions as follows:

${\mathrm{RMSD}}^2={\mathbf{T}}^2+\frac{4}{W}{\mathbf{q}}^T\mathbf{I}\mathbf{q}+2{\mathbf{T}}^T(\mathbf{R}-{\mathbf{E}}_3)\mathbf{C}$$\operatorname{RMSD}^{2}=\mathbf{T}^{2}+\frac{4}{W} \mathbf{q}^{T} \mathbf{I q}+2 \mathbf{T}^{T}\left(\mathbf{R}-\mathbf{E}_{3}\right) \mathbf{C}$

With ${\mathbf{E}}_{\mathbf{3}}$$\mathbf{E_{3}}$ being the $3\times 3$$3 \times 3$ identity matrix, $C$$C$ represents the centre of mass $\frac{1}{W}{\{\sum w_i{x}_i,\sum w_i{y}_i,\sum w_i{z}_i\}}^T$$\frac{1}{W}\left\{\sum w_{i} x_{i}, \sum w_{i} y_{i}, \sum w_{i} z_{i}\right\}^{T}$, and I denotes the inertia tensor.

We should mention that it is practical to work in the center-of-mass reference frame where $C=0$$C=0$. Thus, in this frame, the RMSD can be expressed with fewer arithmetic operations as

${\mathrm{RMSD}}^2={\mathbf{T}}_{\text{COM}}^2+\frac{4}{W}{\mathbf{q}}^T{\mathbf{I}}_{\text{COM}}\mathbf{q}$$\mathrm{RMSD}^{2}=\mathbf{T}_{\text{COM}}^{2}+\frac{4}{W}\mathbf{q}^{T} \mathbf{I}_{\text{COM}}\mathbf{q}$

RMSD for flexible molecules modelled with collective motions

---

We are expanding on our previous work ^[7] by incorporating molecular flexibility through linear collective motions. These can be computed using techniques such as Normal Mode Analysis (NMA) or Principal Component Analysis (PCA). More specifically, let ${\left\{f_i^j\right\}}_N^M$ $\left\{f_{i}^{j}\right\}_{N}^{M}$ be a set of $M$ $M$ vectors that describe the linear collective motions applied to a molecule, where ${\mathbf{f}}_i^j=$ $\mathbf{f}_{i}^{j}=$ ${\{f_{ix}^j,f_{iy}^j,f_{iz}^j\}}^T$ $\left\{f_{i x}^{j}, f_{i y}^{j}, f_{i z}^{j}\right\}^{T}$ . Here, $\mathrm{i}$ $\mathrm{i}$ is the atom index ranging from 1 to $N$ $N$ , and $\mathrm{j}$ $\mathrm{j}$ is the index of the collective motions ranging from 1 to $M$ $M$ . Let ${\left\{{\mu }^j\right\}}^M$ $\left\{\mu^{j}\right\}^{M}$ be the amplitudes of the collective motions for the reference conformation of the molecule. Then the reference flexible coordinates $A^F={\left\{a_i^F\right\}}_N$ $A^{F}=\left\{a_{i}^{F}\right\}_{N}$ are given as ${\mathbf{a}}_i^F={\mathbf{a}}_i+\sum _{j=1}^M{\mu }^j{\mathbf{f}}_i^j$ $\mathbf{a}_{i}^{F}=\mathbf{a}_{i}+\sum_{j=1}^{M} \mu^{j} \mathbf{f}_{i}^{j}$ . Now, to compute the flexible coordinates of the target conformation $A^{\prime F}$ $A^{\prime F}$ , we first add flexible displacements to the rigid coordinates of the reference conformation and then apply the rigid-body transformation to the result. Let ${\left\{{\lambda }^j\right\}}^M$ $\left\{\lambda^{j}\right\}^{M}$ be the amplitudes of the collective motions for the target conformation of the molecule, where the collective motion vectors are the same as in the reference conformation. Similarly to the rigid-body case, let $\mathbf{R}$ $\mathbf{R}$ and $\mathbf{T}$ $\mathbf{T}$ be the rotation matrix and the translation vector of the target rigid-body transformation, respectively, so that the flexible target coordinates ${\left\{{\mathbf{a}}_i^{\prime F}\right\}}_N$ $\left\{\mathbf{a}_{i}^{\prime F}\right\}_{N}$ are given as ${\mathbf{a}}_i^{\prime F}=\mathbf{R}({\mathbf{a}}_i+\sum _{j=1}^M{\lambda }^j{\mathbf{f}}_i^j)+\mathbf{T}$ $\mathbf{a}_{i}^{\prime F}=\mathbf{R}\left(\mathbf{a}_{i}+\sum_{j=1}^{M} \lambda^{j} \mathbf{f}_{i}^{j}\right)+\mathbf{T}$ . Then, the weighted RMSD between positions ${\mathbf{A}}^F$ $\mathbf{A}^{F}$ and ${\mathbf{A}}^{\prime F}$ $\mathbf{A}^{\prime F}$ is given as follows:

${\mathrm{RMSD}}^2(A^F,A^{\prime F})=\frac{1}{W}\sum _i{w}_i{|{\mathbf{a}}_i+\sum _j{\mu }^j{\mathbf{f}}_i^j-\mathbf{R}({\mathbf{a}}_i+\sum _j{\lambda }^j{\mathbf{f}}_i^j)-\mathbf{T}|}^2$ $\operatorname{RMSD}^{2}\left(A^{F}, A^{\prime F}\right)=\frac{1}{W} \sum_{i} w_{i}\left|\mathbf{a}_{i}+\sum_{j} \mu^{j} \mathbf{f}_{i}^{j}-\mathbf{R}\left(\mathbf{a}_{i}+\sum_{j} \lambda^{j} \mathbf{f}_{i}^{j}\right)-\mathbf{T}\right|^{2}$

We can rewrite the previous expression using the quaternion representation of vectors ${\mathbf{a}}_i,\mathbf{T}$ $\mathbf{a}_{i}, \mathbf{T}$ , and the rotation matrix $\mathbf{R}$ $\mathbf{R}$ as

${\mathrm{RMSD}}^2=\frac{1}{W}\sum _i{w}_i{|[0,{\mathbf{a}}_i+\sum _j{\mu }^j{\mathbf{f}}_i^j]-\hat{Q}[0,{\mathbf{a}}_i+\sum _j{\lambda }^j{\mathbf{f}}_i^j]{\hat{Q}}^{-1}-[0,\mathbf{T}]|}^2$ $\operatorname{RMSD}^{2}=\frac{1}{W} \sum_{i} w_{i}\left|\left[0, \mathbf{a}_{i}+\sum_{j} \mu^{j} \mathbf{f}_{i}^{j}\right]-\hat{Q}\left[0, \mathbf{a}_{i}+\sum_{j} \lambda^{j} \mathbf{f}_{i}^{j}\right] \hat{Q}^{-1}-[0, \mathbf{T}]\right|^{2}$

Here, the unit quaternion $\hat{Q}$ $\hat{Q}$ corresponds to the rotation matrix $\mathbf{R}$ $\mathbf{R}$ . Since the norm of a quaternion does not change if we multiply it by a unit quaternion, we may right-multiply the kernel of the previous expression by $\hat{Q}$ $\hat{Q}$ to obtain

${\mathrm{RMSD}}^2=\frac{1}{W}\sum _i{w}_i{|[0,{\mathbf{a}}_i+\sum _j{\mu }^j{\mathbf{f}}_i^j]\hat{Q}-\hat{Q}[0,a_i+\sum _j{\lambda }^j{\mathbf{f}}_i^j]-[0,\mathbf{T}]\hat{Q}|}^2$ $\operatorname{RMSD}^{2}=\frac{1}{W} \sum_{i} w_{i}\left|\left[0, \mathbf{a}_{i}+\sum_{j} \mu^{j} \mathbf{f}_{i}^{j}\right] \hat{Q}-\hat{Q}\left[0, a_{i}+\sum_{j} \lambda^{j} \mathbf{f}_{i}^{j}\right]-[0, \mathbf{T}] \hat{Q}\right|^{2}$

Using the scalar-vector representation of a quaternion, $\hat{Q}=[s,q]$ $\hat{Q}=[s, q]$ , we rewrite the previous RMSD expression as

Performing scalar and vector products in the above equation, then grouping the terms that depend on atomic positions together, and after introducing the inertia tensor $\mathbf{I}$$\mathbf{I}$, the center of mass vector $\mathbf{C}$$\mathbf{C}$ and reintroducing the rotation matrix $\mathbf{R}$$\mathbf{R}$, we obtain

Here, $\mathrm{Tr}()$$\operatorname{Tr}()$ is the matrix trace operator, $D^j$$D^{j}$ is the set of $M3\times 3$$M 3 \times 3$ matrices of cross-products

is the set of $M^{2}3\times 3$$M^{2} 3 \times 3$ matrices of cross-products

and ${\mathbf{B}}^j=\frac{1}{W}{\{\sum w_i{f}_{ix}^j,\sum w_i{f}_{iy}^j,\sum w_i{f}_{iz}^j\}}^T$$\mathbf{B}^{j}=\frac{1}{W}\left\{\sum w_{i} f_{i x}^{j}, \sum w_{i} f_{i y}^{j}, \sum w_{i} f_{i z}^{j}\right\}^{T}$ are the centers of the collective motions. Again, it is practical to choose the reference frame of the target molecule such that $\mathrm{C}=0$$\mathrm{C}=0$. Also, commonly used collective motions, e.g. those computed using the NMA or the PCA, possess the weight orthonormality property, i.e. $\mathrm{Tr}\left({\mathbf{F}}^{jk}\right)={\delta }_{jk}/W$$\operatorname{Tr}\left(\mathbf{F}^{j k}\right)=\delta_{j k} / W$. From now on, we will consider only this type of motions. Thus, in this case, the RMSD equation simplifies to

This is our primary RMSD equation and the main result of this work, we call it master equation. It comprises the rigid contribution ${\mathbf{T}}^2+\frac{4}{W}{\mathbf{q}}^T\mathbf{I}\mathbf{q}$$\mathbf{T}^{2}+\frac{4}{W} \mathbf{q}^{T} \mathbf{I q}$, the flexible contribution $\frac{1}{W}\sum _j({\mu }^{j2}+{\lambda }^{j2})$$\frac{1}{W} \sum_{j}\left(\mu^{j 2}+\lambda^{j 2}\right)$, and the cross-terms. Once the matrices and vectors $\mathbf{I},{\mathbf{B}}^j,{\mathbf{D}}^j$$\mathbf{I},\mathbf{B}^{j},\mathbf{D}^{j}$, and ${\mathbf{F}}^{jk}$$\mathbf{F}^{j k}$ that depend on the number of atoms and the number of collective motions are computed, the calculation time of RMSD is independent of the number of atoms and is at most quadratic with the number of collective motions. In the following, we will explicitly consider several special cases of simplified motions that also simplify the master equation and reduce its computational cost.

RMSD corresponding to a rigid reference conformation

---

A practical implication of the master equation is the expression of RMSD for a flexible target conformation with respect to a rigid reference conformation. In this scenario, all the amplitudes of the collective motions for the reference conformation ${\left\{{\mu }^j\right\}}^M$$\left\{\mu^{j}\right\}^{M}$ are zero, and the RMSD expression simplifies to the following:

In this case, the calculation of one RMSD takes linear time with the number of collective motions $M$$M$.

RMSD corresponding to a pure flexible motion

---

When the study is limited to the flexible movements of a molecule, such as in the refinement of docking poses or the generation of pseudo-random structural ensembles, the master equation becomes simplified

${\mathrm{RMSD}}^2=\frac{1}{W}\sum _j{({\mu }^j-{\lambda }^j)}^2$$\operatorname{RMSD}^{2}=\frac{1}{W} \sum_{j}\left(\mu^{j}-\lambda^{j}\right)^{2}$

Indeed, in this case, the rotation matrix $\mathbf{R}$$\mathbf{R}$ is identity and the translation vector $\mathbf{T}$$\mathbf{T}$ is zero. Thus, the RMSD expression becomes linear in the number of collective motions $M$$M$.

RMSD corresponding to the relative rigid-body motion

---

The RMSD master equation can also be adapted for the particularly useful clustering application case. Here, we compare two possible target conformations, for which the transformation operators are defined with respect to the original reference conformation. Let $A={\left\{a_i\right\}}_N$$A=\left\{a_{i}\right\}_{N}$ define the reference coordinates and $A_1={\left\{a_i^{(1)}\right\}}_N$$A_{1}=\left\{a_{i}^{(1)}\right\}_{N}$ and $A_2={\left\{a_i^{(2)}\right\}}_N$$A_{2}=\left\{a_{i}^{(2)}\right\}_{N}$ be the two target conformations we want to compare. Let set ${\left\{{\lambda }^j\right\}}^M$$\left\{\lambda^{j}\right\}^{M}$ define the collective motion amplitudes of $A_1$$A_{1}$, and set ${\left\{{\mu }^j\right\}}^M$$\left\{\mu^{j}\right\}^{M}$ define the ones of $A_2$$A_{2}$. Let also $R_1$$R_{1}$ and $T_1$$T_{1}$ describe the rigid-body transformation applied to A to obtain $A_1$$A_{1}$, and $R_2$$R_{2}$ and $T_2$$T_{2}$ be the rigid-body transformation applied to $\mathrm{A}$$\mathrm{A}$ to obtain $A_2$$A_{2}$. Finally, let a unit quaternion $[s_{12},q_{12}]$$\left[s_{12}, q_{12}\right]$ correspond to the relative rotation ${\mathbf{R}}_{12}\equiv {\mathbf{R}}_2^T{\mathbf{R}}_1$$\mathbf{R}_{12} \equiv \mathbf{R}_{2}^{T} \mathbf{R}_{1}$ and let the relative translation be ${\mathbf{T}}_{12}\equiv {\mathbf{R}}_1^T({\mathbf{T}}_2-{\mathbf{T}}_1)$$\mathbf{T}_{12} \equiv \mathbf{R}_{1}^{T}\left(\mathbf{T}_{2}-\mathbf{T}_{1}\right)$. Then, the weighted RMSD between $A_1$$A_{1}$ and $A_2$$A_{2}$. is given by a generalized version of the master equation as follows:

Computational analysis revealed that the Root Mean Square Deviation (RMSD) of the CTP molecule in relation to the Ncds-2 protein is 2.68, considering 43 to 43 atoms. In contrast, the RMSD of the NMN molecule with respect to the Ncds-2 protein is 2.32, considering 30 to 30 atoms.

Conclusion

---

Molecular docking simulation is a computational process that models the interaction between two or more molecules. This simulation aids in understanding the binding mechanisms of molecules and the strength of these interactions. Through molecular docking simulation, we have computationally validated the catalytic activity of Ncds towards the synthesis of $\mathrm{NCD}$$\mathrm{NCD}$, thereby providing valuable theoretical guidance for subsequent wet-lab experiments.

Reference

---

[1] Raj, A., & van Oudenaarden, A. (2008). Nature, Nurture, or Chance: Stochastic Gene Expression and Its Consequences. Cell, 135(2), 216–226.

[2] Guo, X., Liu, Y., Wang, Q., Wang, X., Li, Q. X., Liu, W-J., & Zhao, Z. K. (2020). Non-natural Cofactor and Formate‐Driven Reductive Carboxylation of Pyruvate. Angewandte Chemie International Edition, 59(8), 3143–3146.

[3] Ji, D., Wang, L., Hou, S., Liu, W., Wang, J., Wang, Q., & Zhao, Z. K. (2011). Creation of Bioorthogonal Redox Systems Depending on Nicotinamide Flucytosine Dinucleotide. Journal of the American Chemical Society, 133(51), 20857–20862.

[4] Wang, J., Guo, X., Wan, L., Liu, Y., Xue, H., & Zhao, Z. K. (2022). Engineering Formaldehyde Dehydrogenase from Pseudomonas putida to Favor Nicotinamide Cytosine Dinucleotide. Chembiochem: A European Journal of Chemical Biology, 23(7), e202100697.

[5] Zelle, RM., de Hulster, E., van Winden, WA., de Waard, P., Dijkema, C., Winkler, AA., Geertman, J-MA., van Dijken, JP., Pronk, JT., & van Maris, AJA. (2008). Malic Acid Production by Saccharomyces cerevisiae: Engineering of Pyruvate Carboxylation, Oxaloacetate Reduction, and Malate Export. Applied and Environmental Microbiology, 74(9), 2766–2777.

[6] Borda-Molina, D., José Salvador Montaña, María Mercedes Zambrano, & Baena, S. (2017). Mining lipolytic enzymes in community DNA from high Andean soils using a targeted approach. Antonie Van Leeuwenhoek International Journal of General and Molecular Microbiology, 110(8), 1035–1051.

[7] Neveu, E., Popov, P., Hoffmann, A., Migliosi, A., Besseron, X., Grégoire Danoy, Pascal Bouvry, & Sergei Grudinin. (2018). RapidRMSD: rapid determination of RMSDs corresponding to motions of flexible molecules. Bioinformatics, 34(16), 2757–2765.