Predicting Protein Folding Pathways With Quadratic Constraints on Rates of Entropy Change: A Nonlinear Optimization-Based Control Approach Open Access

The protein folding problem can be understood in two key parts [1,2], which build upon Anfinsen's seminal work [3] and Levinthal's proposition [4]. In the former, the aim is to predict a protein's three-dimensional (3D) structure solely from its amino acid sequence. In the latter, the objective is to compute transient conformations of proteins during folding (see Fig. 1).

Understanding protein folding as a guided process involving specific pathways or mechanisms stems from Levinthal's paradox [4]. According to this paradox, protein folding through a random search of every possible protein conformation is implausible. This is because the number of potential protein configurations is astronomically vast and a random search would lead to an incredibly long folding time (Levinthal's time). Despite the inherent challenge of Levinthal's paradox, incorporating a minimal energy bias against energetically unfavorable protein conformations can significantly reduce the folding time scale to biologically relevant regimes [5]. Furthermore, the difference in entropy between denatured and folded proteins contributes to the complexity of folding. In particular, denatured proteins exhibit a markedly higher conformational entropy compared to their folded native state. Moreover, two-state proteins preferentially adopt pathways that minimize entropy loss during folding [6].

The limited temporal resolution of protein structural data presents a significant challenge to experimentally elucidating protein folding pathways [7]. Consequently, algorithmic prediction of protein folding pathways has emerged as a powerful complementary approach, finding applications in computer-aided drug discovery [8] and studies of protein misfolding diseases [9]. Furthermore, folding simulations hold promise for engineering applications, enabling quantitative prediction of mobility and range of motion in protein-based nanomachines [10,11].

One category of widely used algorithms for predicting the final folds of proteins is provided by what is known as knowledge-based approaches. These algorithms, which are based on pattern recognition and machine learning, infer the 3D structures of folded protein shapes by harnessing the predictive power of the amino acid sequences (akin to Anfinsen's dogma) and leveraging extensive datasets of existing folds [12]. For instance, AlphaFold's core component is a neural network trained on the protein data bank (PDB) structures. This network predicts probability distributions of inter-residue atom distances and leverages the patterns derived from the vast protein sequence landscape through the integration of sequence databases and multiple sequence alignments [13]. However, knowledge-based approaches cannot describe the interactions between proteins and nucleic acids. Furthermore, they fail to take into account the stability of the protein molecule and the rates at which the folding process occurs, i.e., the folding kinetics (see, e.g., Ref. [14]).

In contrast to knowledge-based strategies, physics-based methods utilize basic physical principles to simulate the process of folding while computing the transient and final 3D structures of protein molecules [15]. While physics-based strategies supply trustworthy data concerning the temporary structures in the folding process, these simulations require significant processing time and specialized hardware (e.g., Anton family of supercomputers for atomic-level folding simulations [16]).

To overcome the computational burden of physics-based protein folding algorithms, the innovative kinetostatic compliance method (KCM) framework has been developed in the literature [17–20]. KCM models protein molecules as a multitude of rigid nanolinkages forming a nanokinematic mechanism subject to movement limitations imposed by chemical bonds. To address the significant computational expense associated with all-atom molecular dynamics, KCM iteratively refines protein structures by varying their dihedral angles in a kinetostatic manner. These angles, which define the relative orientations of bonded atoms, are varied under the influence of nonlinear interatomic force fields through an iterative process driving the protein molecule conformation toward a minimized energy state corresponding to the final folded structure.

Kinetostatic compliance method has been a valuable tool for understanding hydrogen bond formation and its influence on protein mobility since its introduction [21]. Moreover, KCM has facilitated the design of peptide-based nanomachines [11]. However, a critical knowledge gap exists: the conventional KCM does not consider how it perturbs the rate of entropy change in kinetostatic protein folding.

The idea of constraining protein folding dynamics using a guiding control input finds its roots in targeted molecular dynamics (TMD) [22,23]. In TMD, a proportional control input regulates the distance between a protein molecule and its target conformation, aiming to expedite the transition to the final folded state.

One of the few attempts to constrain the folding simulations by synthesizing proper guiding control inputs is due to Arkun and collaborators [24,25]. They utilized a bead-spring model of the protein molecule, assuming a linear time-invariant system. They subsequently employed a linear quadratic regulator (LQR) control scheme to guide the folding process. Their controller aimed to navigate away from high-energy regions of the folding landscape while simultaneously regulating the rate of molecule entropy change. However, the simplistic force fields used by Refs. [24] and [25] cannot capture the nonlinear and complicated dynamics of proteins under the electrostatic and van der Waals interatomic effects.

To address the limitations of the protein bead-spring LQR-based folding control framework due to Arkun and collaborators [24,25], we recently proposed a control-theoretic interpretation of KCM [26]. We demonstrated that constraints on the rate of change of entropy can be incorporated into the KCM folding dynamics through control inputs that steer the folding simulation. These control inputs are obtained by solving a large-scale quadratic program (QP) with constraints on a physically relevant norm of the control input. While our initial work utilized interior point algorithms to solve box-constrained QPs for computing the folding control inputs, it encountered a significant limitation. Specifically, it relied on polytopic inequalities, as noted in Ref. [26], which are compatible with rapid numerical solvers (see, e.g., Ref. [27]). This approach contrasts with more realistic quadratic constraints on the inputs. Notably, standard methods in biochemistry, which represent protein configurational entropy, employ quadratic representations based on ellipsoidal norms [28].

This paper makes two key contributions to the KCM-based protein folding framework. First, we address a significant gap in KCM-based protein folding literature. While our previous work constrained the rates of entropy change in the KCM paradigm for the first time [26], it suffers the shortcoming of encoding these constraints via polytopic inequalities with an inevitable loss of accuracy in predicting the folding pathways. To overcome this challenge, we propose a novel approach that employs a computationally efficient trust-region subproblem (TRS) algorithm for constraining the rates of entropy change via ellipsoidal inequalities (Algorithm 1). Moreover, this algorithm offers a significant advantage by eliminating the need for outer iterations, thereby enhancing computational efficiency. This computationally efficient control synthesis algorithm, which is taken from the TRS optimization literature [29], is based on solving a generalized eigenvalue problem (GEP). This approach ensures that the TRS-based folding control inputs induce a closed-loop dynamics faithful to the KCM reference vector field while satisfying the ellipsoidal constraints on the folding control input.

Our second contribution focuses on leveraging the large-scale GEP-based TRS solver developed by Adachi et al. [29] to generate optimization-based control inputs. While Adachi et al. primarily targeted the optimization literature, where solution continuity was not a primary concern [29], control systems considerations necessitate investigating continuity properties of the generated control inputs (see, e.g., Ref. [30] for humanoid robots). Remarkably, the application of Adachi et al.'s solver and its continuity properties within the context of nonlinear optimization-based control remains unexplored. We address this gap by investigating the piecewise continuity of pointwise optimal control laws obtained from the GEP-based TRS solver. We also establish conditions for existence and uniqueness of closed-loop trajectories under these control inputs (Propositions 5.2 and 5.4). Other contributions of this paper include the novel derivation of KCM protein folding dynamics from a Lagrange–Rayleigh formalism (Sec. 3.1) and the first formal proofs of convergence for the KCM folding scheme to folded conformations in the literature (Sec. 3.2).

This paper significantly expands upon a preliminary version presented as a conference poster [31]. In this paper, we delve deeper into the control synthesis algorithm by providing a detailed description, discussing its computational complexity, and exploring its numerical properties. Additionally, we present a comprehensive analysis (absent from the conference poster) investigating the piecewise continuity of control inputs and establishing conditions for existence and uniqueness of closed-loop protein folding trajectories. Finally, the algorithm effectiveness is rigorously assessed through more extensive numerical simulations.

The remainder of this paper is structured as follows. First, in Sec. 2, we briefly present the nanokinematic structure of protein molecules and the main underpinnings of the KCM framework for simulating the kinetostatic folding of proteins. We further demonstrate how the KCM-based approach to protein folding can be formulated as a control synthesis problem. Thereafter, we formulate the main control problem in Sec. 3 and present a pointwise optimal control scheme for solving it in Sec. 4. We subsequently present conditions for existence and uniqueness of solutions to the resulting closed-loop protein folding dynamics in Sec. 5. After presenting the simulation results in Sec. 6, we conclude the paper with final remarks, the impact of the current work beyond protein folding, and future directions in Sec. 7.

Throughout this paper, $R+$ and $Z0+$ denote the set of all non-negative real numbers and non-negative integers, respectively. The set $S1$ denotes the unit circle. The set $SO(3)$ denotes the special orthogonal group of order three. The matrix $In∈Rn×n$ denotes the identity matrix. Given the vector $x∈Rn$ and symmetric positive definite matrix $Y∈Rn×n$⁠, $|x|2:=x⊤x$ denotes the vector Euclidean norm, $|x|∞:=maxi|xi|$ denotes the vector $∞$-norm, and $|x|Y:=x⊤Yx$ denotes the vector ellipsoidal norm. Given the integer n, $In$ denotes the identity matrix of size n. Hence, $|x|In=|x|2$⁠. Given a matrix $X$⁠, $R(X)$ and $N(X)$ denote the range and null space of the matrix, respectively. Additionally, $X≻0$ (⁠$X≽0$⁠) indicate that $X$ is a positive definite (resp., semidefinite) matrix. Given a set $S⊆Rn$⁠, $co(S)$ (⁠$co ¯(S)$⁠) denote the convex hull (respectively, convex closure) of the set. Given a smooth map $h(·)$ from an open subset of $Rm$ to an open subset of $Rn$⁠, $dhx$ denotes the Jacobian of $h(·)$ at $x$⁠. If $h(·)$ is a real-valued function, $dhx⊤=∇xh(x)$⁠, where $∇xh(x)$ denotes the gradient of $h(·)$ at $x$⁠. Given a vector field $Θ:Rn→Rn$⁠, $FI[Θ](·)$ is the Filippov set-valued map associated with the vector field.

In this section, we provide a summary of the KCM framework for modeling the kinetostatic folding of proteins in vacuo.

Proteins are macromolecules with complex structures and intricate dynamics governing their folding process. These long molecular chains consist of peptide planes that are interconnected by chemical bonds. The functionality of protein molecules depends critically on their 3D structure (i.e., their conformation), which can be directly deduced from the linear sequence of amino acids forming their polypeptide chain.

A protein backbone chain schematic is depicted in Fig. 2, where each individual peptide plane is formed by six coplanar atoms that are covalently bonded (visualized as red lines in Fig. 2). The coplanarity assumption about the $α$-carbon $Cα$⁠, carbonyl CO, amide nitrogen NH, and another $α$-carbon is motivated by the wealth of high-resolution X-ray crystallographic data [2]. This fundamental assumption underpins the representation of protein molecules as intricate nanomechanisms with numerous degrees-of-freedom (DOFs) (see, e.g., Refs. [18], [26], and [32–34]). Accordingly, the key building blocks of the protein kinematic chain are the peptide planes that function as nanoscale rigid linkages.

The rotation of the nanolinkages in the protein kinematic chain with respect to each other is facilitated by the alpha-Carbon atoms (Fig. 2), which play the role of hinges connecting peptide planes together. The first and the last $Cα$ atoms are bonded to an N-terminus–peptide plane and a C-terminus–peptide plane, respectively.

Represents the configuration of the backbone chain. The dihedral angle vector $θ$ resides within the 2N-dimensional configuration space $Q=S1×⋯×S1$⁠, where each individual unit circle $S1$ corresponds to a single dihedral angle within the backbone chain.

Each dihedral angle in the conformation vector $θ$ in Eq. (1) corresponds to a single DOF in the protein molecule's kinematic chain. For each DOF, we can define a unit vector, $uj$⁠, $1≤j≤2N$⁠. These vectors align with the axes of rotation around which the nanokinematic chain can rotate. As illustrated in Fig. 2, the unit vectors $u2i$ and $u2i+1$ represent the directions along two crucial bonds within the ith peptide plane: the $Cα−C$ bond and the $N−Cα$ bond. These bonds define the primary rotation axes within the chain's segments. Serving as key structural references, $u1$ and $u2N$ denote the unit vectors associated with the amino terminus (N-terminus) and carboxy terminus (C-terminus), respectively. The reference conformations are established by setting all dihedral angles to zero, i.e., $θi=0$ for $1≤i≤2N$⁠. This configuration results in the protein adopting a linear coil conformation, representing a denatured state. In this state, the protein is stretched or completely unfolded, lacking the specific three-dimensional structure that is characteristic of its functional, native form. Such a state is typically induced by various physical or chemical means, including heat, changes in pH, Chaotropic agents like urea, or mechanical stretching.

While dihedral angles define the rotational DOFs within a protein's kinematic chain, additional information is necessary to fully capture the spatial orientation of the rigid peptide nanolinkages. This information comes in the form of body vectors, denoted by $bj$⁠, $1≤j≤2N$⁠. These vectors comprehensively define the relative spatial arrangement of coplanar atoms within each peptide plane. Consequently, any relative positional relationship between two atoms within a plane can be mathematically expressed as a linear combination of their corresponding body vectors, i.e., $k1mb2i+k2mb2i+1$⁠, where the coefficients $k1m$ and $k2m$⁠, $1≤m≤4$⁠, are constants specific to the atom pair and are identical throughout the peptide chain (see, e.g., Refs. [18] and [19] for more details). The body vector $b2i$ is oriented along the line connecting the alpha Carbon atom to the Nitrogen atom on the ith peptide plane, while $b2i+1$ aligns with the chemical bond linking the Nitrogen atom of the ith peptide plane to the alpha Carbon of the following peptide plane. Due to their distinct directional alignments within the peptide plane structure, $b2i$ and $b2i+1$ are inherently linearly independent (Fig. 2). The body vectors $bj$ and unit vectors $uj$ act in concert to provide a comprehensive depiction of the protein molecule's conformation as it dynamically varies with respect to the vector $θ$⁠, which encompasses the peptide dihedral angles.

Kinetostatic compliance method framework leverages the established experimental observation that it is possible to capture the essence of protein folding dynamics while neglecting inertial forces (see, e.g., Refs. [24], [25], and [35]). According to KCM, protein dihedral angles evolve kinetostatically under the influence of interatomic force fields.

Consider a peptide chain composed of $Na$ atoms and $N−1$ peptide planes. The dihedral angle vector, denoted by $θ$ defined in Eq. (1), encodes the conformational state of the chain. The Cartesian coordinates of any two atoms $ai$⁠, $aj$ within the chain are represented by $ri(θ)$⁠, $rj(θ)$⁠, respectively (Eq. (4)). Their Euclidean distance is then calculated as $dij(θ):=|ri(θ)−rj(θ)|2$⁠. The specific parameters pertaining to atom charges, van der Waals radii, interatomic distances, dielectric constant, potential well depths, and force weights are readily available in Ref. [19] and the provided references.

where $Gelec(θ)$ and $Gvdw(θ)$ represent the protein's electrostatic potential energy and van der Waals interatomic potential energy, respectively (detailed expressions can be found in Ref. [18]). Consequently, the net forces acting on each atom $ai$⁠, $1≤i≤Na$⁠, due to both Coulombic and van der Waals interactions, can be obtained from the negative gradients of the respective individual energy terms: $Fielec(θ)=−∇riGelec$ and $Fivdw(θ)=−∇riGvdw$⁠. It is remarked that the molecular dynamics formulations begin with Cartesian coordinates of the atoms given by $ri$ for $1≤i≤Na$⁠, before transitioning to dynamics described in generalized coordinates (see, e.g., Ref. [36]). In the KCM framework, the generalized coordinates are given by the dihedral angles $θ1,…,θ2N$⁠. Hence, $ri=ri(θ1,…,θ2N)$ for $1≤i≤Na$⁠. The free energy of the molecule $G(θ1,…,θ2N)$ is indeed a function of $θ$⁠, and thus, should explicitly relate to the Cartesian coordinates as $G(r1(θ1,…,θ2N),…,rNa(θ1,…,θ2N))$⁠.

where $J(θ)∈R6N×2N$ represents the molecule chain Jacobian at conformation $θ$⁠. This matrix translates the 6N-dimensional force vector into the corresponding 2N-dimensional torque vector governing dihedral angle evolution (see Ref. [18] for detailed derivations). Consequently, the changes in dihedral angles at each protein conformation take place under the direct influence of the kinetostatic folding torque vector acting on the peptide backbone. Specifically, the vector $τ(θ)∈R2N$⁠, which aligns with the steepest-descent direction of the total free energy $G(θ)$ in the protein conformation landscape, guides the dihedral angle variations during folding.

Within the context of protein folding, at each local minimum $θ∗$ of the aggregate free energy function $G(θ)$ defined in Eq. (5), the corresponding torque vector $τ(θ∗)$ vanishes identically. This signifies the protein molecule's kinetostatic stationarity at folded conformations, where the absence of net kinetostatic torque reflects the balanced internal forces within the protein structure.

where $τ(θk)=J⊤(θk)F(θk)$ due to Eq. (6). The positive constant $κ0$ is chosen sufficiently small to prevent excessive angle variations and is heuristically optimized.

A flowchart of the successive kinetostatic folding iteration is depicted in Fig. 3. In the dihedral angle update iterations given by (7), the normalized torque vector $τ(θk)|τ(θk)|∞$ governs the incremental updates to the dihedral angles at each conformation $θk$ (see Ref. [20] for a sign gradient descent-based version of this iteration). The iterative process governed by Eq. (7) continues until the molecule's aggregated free energy $G(θ)$ converges to the vicinity of a local minimum within the free-energy landscape. Convergence is achieved numerically when the kinetostatic folding torque vector norm falls below a predefined tolerance $τtol>0$⁠, i.e., when $|τ(θk∗)|2<τtol$⁠, for some positive integer $k∗$⁠. The protein folding pathway transient conformations are given by ${θk}k=0k∗$⁠.

In this section we state the smoothness properties of the aggregated energy of the protein molecule, $G(θ)$⁠, in the vicinity of folded conformations (see, e.g., Ref. [37] for further details). In particular, in a vicinity of a folded conformation $θ=θ∗$⁠, $G(θ)$ is a smooth function of the dihedral angle vector $θ$ and assumes a local minimum in this folded conformation. Therefore, $θ∗$ is a critical point of $G(·)$ and $∇θG(θ∗)=0$⁠. Finally, in the so-called protein metastable states including the folded conformations, the Hessian of $G(·)$ is positive definite.

links the gradient of the aggregated free energy to the generalized torques acting on dihedral angles, as demonstrated by Noguti and Gō [38]. Equation (8) is applicable to noncyclic molecules including proteins.

In this section, we recast the KCM-based folding numerical scheme presented in Sec. 2.2 as a nonlinear affine control system. Our point of departure is the Lagrange–Rayleigh formalism resulting in an Euler–Lagrange control system. This controls-oriented framing, which first appeared in our earlier work [26,34], draws inspiration from the field of targeted molecular dynamics (TMD) [22,23] and Arkun et al.'s optimal control framework, where they employed linear spring-like forces described by a connectivity matrix [24,25,39]. In contrast to Arkun et al., our KCM-based line of work utilizes a nonlinear dynamical system to capture the dynamics within the molecule's dihedral angle space.

where the Jacobian matrix, $J(θ)$⁠, and the generalized force vector, $F(θ)$⁠, are identical to those defined in Eq. (6). Equation (10) describes the dynamics of the dihedral angles of a protein molecule subject to interatomic forces and a folding control input, $uc$⁠. In the KCM dynamics, the folding control input, $uc$⁠, can represent either physical forces, such as those exerted by optical tweezers or solvation effects, or nonphysical fictitious forces designed to enhance the convergence rate toward folded conformations in simulations performed in vacuo. In the case of the latter, which is the subject of this paper, the folding control input functions similarly to the forces used in targeted molecular dynamics.

Notably, employing a forward Euler scheme with an integration step equal to $κ0$ recovers the KCM iteration in Eq. (7) from the reference vector field in Eq. (12). In summary, the successive kinetostatic folding iteration in Eq. (7) is related to the reference vector field given by Eq. (12). This vector field is induced by applying the folding control input (11) to the KCM-based nonlinear control affine system in Eq. (10). Furthermore, the trajectories of the vector field in Eq. (12) provide the protein folding pathways.

Remark 1. While the controls-oriented perspective in Ref. [26] was not explicitly addressed in the original KCM-based protein folding works (e.g., Refs. [18] and [19]), it offers a valuable foundation for extending the KCM framework. As discussed later, this control-theoretic interpretation allows us to constrain the molecule entropy changes introduced by the folding control input, leading to more realistic pathway predictions.

Remark 2. The KCM reference vector field described by Eq. (12), originally developed through the pioneering efforts of Kazerounian, Ilieş, and their collaborators (e.g., Refs. [18] and [41]), has been widely utilized in the KC esponds to a veM literature. This vector field has demonstrated remarkable success in numerical protein folding simulations and has been effectively utilized in significant engineering applications, such as design of peptide-based nanomechanisms [42]. Unfortunately, its use has lacked formal theoretical guarantees and proofs of convergence to folded conformations. In this paper, we address this shortcoming by considering a very similar reference vector field (see Sec. 3.2) for which convergence and smoothness can be easily established. Our simulations will include comparisons of both reference vector fields.

In an open neighborhood of a specific folded conformation—excluding the conformation itself—the new KCM vector field given by Eq. (13) is smooth, thanks to the smoothness properties of the aggregated free energy of protein molecules (Sec. 3.2). The gradient flow induced by the reference vector field (13) has well-defined absolutely continuous solutions. Furthermore, using the results in Cortés' work [43, Proposition 7 and Theorem 8], it can be shown that the KCM flow $θ˙=ΘKCM2(θ)$ with initial conformations in a neighborhood of a given folded conformation is guaranteed to converge to it in finite time due to positive definiteness of the Hessian of the aggregated free energy of the molecule (Sec. 3.2).

Figures 4 and 5 compare the folding pathways of a protein backbone chain with $N=5$ peptide planes, starting near a $β$-pleated sheet conformation (with simulation numerical parameters detailed in the $β$-pleated sheet simulations in Sec. 6), under two different reference vector fields. The KCM reference vector fields $ΘKCM(θ)$ in Eq. (12) and $ΘKCM2(θ)$ in Eq. (13) are derived from the control inputs $uKCM(θ)$ and $uKCM2(θ)$⁠, as defined in Eqs. (11) and (14), respectively, within the KCM dynamics (10). As illustrated in the figures, although the original reference vector field commonly used in KCM folding studies converges more quickly, it does not possess the continuity properties of the modified vector field. Indeed, the conventional KCM control input starts oscillating after around 6000 iterations. Therefore, we continue our theoretical derivations using the modified reference vector field as outlined in Remark 2 and present the numerical simulation results for both vector fields in the sequel.

In this section we consider the controls-oriented protein folding KCM dynamics and investigate how excessive entropy changes due to the folding control input can be constrained.

We first derive an expression for the change of the molecule entropy during folding. This derivation is motivated by a key observation in the field of structural biochemistry, which states that excessively rapid changes of entropy should be avoided during the guided folding process [6,24,39]. Multiple established methodologies exist to calculate the entropy changes experienced by a protein molecule during folding. These methodologies can be broadly classified into two categories: those based on configurational entropy computations (e.g., Ref. [44]) and those rooted in thermodynamics-derived operational definitions (e.g., Refs. [45] and [46]). In this work, we will employ the Clausius thermodynamic expression outlined in Refs. [45] and [46]. This approach shares similarities with the derivations presented by Arkun et al. [24], as both methods leverage the principles of classical thermodynamics.

where for some symmetric positive definite matrix $B$ and constant $Ω0$⁠. When $B=I2N$⁠, we arrive at $|uc|2≤Ω0$⁠.

As it can be seen from Eq. (20), we establish an ellipsoidal constraint on the control input $uc$⁠. This constraint directly limits the total rate of entropy change $S˙$⁠, rather than just the ‘extra’ or ‘deviating’ component $∇θG(θ)⊤uc$⁠. This approach is justified by Eq. (16), where the ‘natural’ torque term $J⊤(θ)F(θ)$ is influenced by the KCM free-energy gradients. By bounding $|uc|B$⁠, we ensure that $S˙$ remains within a permissible range. Consequently, controlling the magnitude of $uc$ not only restricts $∇θG(θ)⊤uc$ but also imposes a direct limit on $S˙$⁠.

Remark 3. As demonstrated in Eq. (17), the gradient of the aggregated free energy, $∇θG(θ)$⁠, is crucial for determining the physics-based rate of entropy change. Additionally, the control input $uc$ is employed to expedite convergence to folded conformations and to model effects such as solvation or optical tweezer forces, though these are not the focus of this study. The magnitude of this control input can artificially alter the rate of entropy change. To mitigate this, our folding control input design aims to closely align the reference vector field with $ΘKCM2(·)$ in Eq. (13), induced by the control input $uKCM2(θ)$ in Eq. (14), while adhering to norm constraints to prevent significant deviations in the natural rate of entropy change.

We are now ready to state the main control problem addressed in this paper.

Consider the protein folding nonlinear control affine system given by Eq. (10). Furthermore, consider the KCM reference vector field $ΘKCM2(·)$ in Eq. (13) induced by the folding control input in Eq. (14). The Protein Folding Control Problem Subject to Entropy Loss Constraints (PFCP) is concerned with two objectives: (1) Synthesize a folding control input $uc(·)$ such that the resulting closed-loop dynamics is closely aligned with the KCM reference vector field while adhering to the ellipsoidal constraint in Eq. (20); and (2) Find conditions for existence and uniqueness of the trajectories of the closed-loopy dynamics, which are the folding pathways subject to the constraint given by Eq. (20).

To achieve the first objective of the PFCP formulated in Sec. 3.3, we first utilize a so-called optimal decision strategy (ODS) to derive a pointwise optimal control law (see, e.g., Refs. [47] and [48] for the pioneering work on ODS), which happens to be the solution to a TRS. We then present an efficient numerical algorithm for computing this optimal control law at each protein conformation.

The proposed control scheme in this section for solving the PFCP leverages an ODS framework to achieve pointwise optimality.

The free vectors in $S(θ(t))$ give rise to the set of admissible dihedral angle folding rate vectors, denoted by $St(θ(t))$⁠. Specifically, every element within $St(θ(t))$ corresponds to a vector in $S(θ(t))$ after its starting point is translated by $θ(t)$⁠. Hence, under any input $uc(t)$ satisfying the constraints, the feasible tangent vectors along the folding trajectory reside within $St(θ(t))$ (Fig. 6).

The following lemma provides the necessary and sufficient conditions for existence of solutions to Eq. (23).

The pair$(λ∗(θ),uc∗(θ))$satisfying the conditions of Proposition 4.1 is called the TRS solution at$θ$. We provide an efficient GEP-based algorithm for computing this solution in Sec. 4.2.

Remark 4. In the limit where the constraints in Eq. (23) become sufficiently relaxed, or equivalently, as $Ω0→∞$⁠, the generated control inputs converge to the conventional KCM control torque vector $uKCM2(·)$ given by Eq. (14).

To solve the TRS in Eq. (23) for folding control input synthesis, we employ the approach due to Adachi et al. [29], which relies on constructing a generalized eigenvalue problem (GEP).

and denote its largest eigenvalue by $λmax∗$⁠.

Algorithm 1 taken from Ref. [29] solves the TRS in Eq. (23). The algorithm can be explained as follows. The solutions to the TRS in Eq. (23) at each protein conformation $θ$ can be separated into two types: (i) interior solutions with $|uc∗(θ)|B<Ω0$⁠; and, (ii) boundary solutions with $|uc∗(θ)|B=Ω0$⁠.

If we have an interior solution at a protein conformation $θ$⁠, then $uc∗(θ)$ can be obtained by solving the linear system $Quc∗(θ)=−G(θ)$ because $Q≻0$⁠. Accordingly, given a protein molecule conformation $θ$⁠, we first compute the solution candidate $p∗=Q−1G(θ)$⁠. If it satisfies $|p∗|B<Ω0$⁠, then it is indeed the solution to the TRS in Eq. (23) at conformation $θ$⁠. Additionally, due to Eq. (25c) in Lemma 4.1, $λ∗(θ)=0$⁠.

On the other hand, if $|p∗|B≥Ω0$⁠, then we are dealing with a boundary case at $θ$⁠. In this case, one needs to solve for the rightmost eigenvalue of the matrix pencil $M˜θ(λ)$⁠, i.e., $λ∗(θ)$⁠, and its associated eigenvector.

Remark 5 (Complexity of Algorithm 1). The computational cost of Algorithm 1 for computing the interior solutions is that of solving a linear system (e.g., the standard method of using the LU decomposition of $Q$⁠), which is $O(N3)$ in the worst case [49]. As demonstrated by Ref. [29], the computational complexity of the algorithm for boundary solution calculations scales as $O(N3)$ floating-point operations (flops) in the case of dense matrices $Q$ and $B$⁠. However, for large and sparse matrices $Q$ and $B$⁠, the cost becomes effectively equivalent to that of an Arnoldi-type method for computing a single rightmost eigenpair. This latter cost is typically proportional to a constant factor multiplied by the cost of a single matrix–vector product.

Finally, we summarize some of the key findings from Ref. [29] regarding the matrix pencils (26) and (27) in the next lemma. We will leverage these results to prove the important properties of our proposed TRS-based control law in Sec. 5.

Lemma 4.2 [29]. Consider a given$θ∈Q$and assume that$|Q−1G(θ)|B≥Ω0$. Consider the matrix pencils$M˜θ(λ)$and$P(λ)$in Eqs. (26) and (27), respectively. The following properties hold.

1. $M˜θ(λ)$is a regular matrix pencil, namely, the determinant of$M˜θ(·)$is non-zero for some$λ$and the number of its eigenvalues is equal to its size.

2. The eigenvalues of$M˜θ(λ)$are the roots of the rational equation

   $hθ(λ):=∑i=12Nfi(θ)(λ−μi)2−Ω02$

   (28)

   where$fi(·)$are smooth functions of$θ$, and$μi$⁠, $1≤i≤2N$, are the eigenvalues of the matrix pencil$P(λ)$. Furthermore,$hθ(λ)$is strictly decreasing on$(μ2N,∞)$, where$μ2N=λmax∗$is the largest eigenvalue of the matrix pencil$P(λ)$⁠.

3. The largest real eigenvalue of$M˜θ(λ)$is$λ∗(θ)$, where$(λ∗(θ),uc∗(θ))$is the solution pair to the TRS in Eq. (23).

4. $λ∗(θ)>λmax∗$. Additionally,$|y1|B≠0$and$G⊤(θ)y2≠0$inAlgorithm 1.

In this section, we address the second objective of PFCP formulated in Sec. 3.3. Specifically, we investigate the properties of the TRS-based folding control input and provide sufficient conditions for existence and uniqueness of the folding trajectories under this input.

To facilitate subsequent analysis within this section, we leverage the closed-form solution provided by the following proposition. While Algorithm 1 offers a numerically stable scheme for computing the TRS-based pointwise optimal control law in Eq. (23), the closed-form solution is more suitable for further analysis.

We have the following proposition regarding the continuity of the TRS-based pointwise optimal control law given by Eq. (31) on $Eθ+$ and $Eθ−$⁠.

Proposition 5.2. The control input given by Eq. (31) is continuous on the open set$Eθ+∪Eθ−$for all protein conformations$θ∈Q$. Moreover, this control input admits a continuous extension to the closure$E¯θ−=Eθ−∪Eθ$⁠.

which, due to Proposition 5.2, is a piecewise smooth vector field with potential discontinuities on the surface $Eθ$⁠. In particular, whenever $λ∗(θ)≠0$ for some $θ∈Eθ$⁠, the control input continuity will be lost.

To investigate the existence and uniqueness of protein folding trajectories induced by the piecewise continuous vector field (32), we need to invoke the notion of Filippov solutions [50]. In particular, we need to compute the Filippov set-valued map of the vector field in (32) and investigate the associated differential inclusion solutions.

where $FI[ΘODS](·)$ is given by Lemma 5.3. If this solution exists, then it is an absolutely continuous trajectory $θ(·)$ satisfying (34) for almost all time instances in its maximal interval of existence.