Programmable Potentials: Approximate N-body potentials from coarse-level logic

This paper gives a systematic method for constructing an N-body potential, approximating the true potential, that accurately captures meso-scale behavior of the chemical or biological system using pairwise potentials coming from experimental data or ab initio methods. The meso-scale behavior is translated into logic rules for the dynamics. Each pairwise potential has an associated logic function that is constructed using the logic rules, a class of elementary logic functions, and AND, OR, and NOT gates. The effect of each logic function is to turn its associated potential on and off. The N-body potential is constructed as linear combination of the pairwise potentials, where the “coefficients” of the potentials are smoothed versions of the associated logic functions. These potentials allow a potentially low-dimensional description of complex processes while still accurately capturing the relevant physics at the meso-scale. We present the proposed formalism to construct coarse-grained potential models for three examples: an inhibitor molecular system, bond breaking in chemical reactions, and DNA transcription from biology. The method can potentially be used in reverse for design of molecular processes by specifying properties of molecules that can carry them out.


Introduction
Multi-body interactions are ubiquitous in nature and happen at all scales from atomic (quantum description) to molecular (classical approach) to macro scales.A systematic analysis these interactions may unfold the fundamental principles governing a given system.For example, understanding the biophysics of protein folding gives insight into disease pathologies 1 .This understanding can be leveraged to develop new vaccines and drug therapies.Engineering these new products requires accurate and computationally tractable models.
Systems having multibody interactions, in fundamental physics, are often formulated as a "N-body potential" problem.In order to fully understand these systems a large number of experiments are needed.Conducting experiments may be expensive and at times even impossible.Another approach is to analyze the N-body potential governing the system dynamics.However, at the quantum level, it may be difficult to determine these potentials from first principles due to the complexity of the system.The computational complexity for ab initio methods can scale exponentially in the number of electrons, limiting the practical size of the system to a few thousand atoms 2-4 .Even if the detailed potential is determined, it may not be immediately useful.Such is the case when the properties or behaviors of interest are at a coarser level than that of the detailed potential and simulating the detailed dynamics is too expensive.Very coarse approaches such as those of master equation 5 lack predictability on molecular spatial and time scales due to the assumptions with which they are derived.A potential that models the system is required if one is to make predictions about the system.
It is profitable to restrict one's efforts to considering approximate potentials that respect known behavior.Such coarse-level descriptions may be determined from experimental observation and may correspond to trajectories in some transformed (reaction) coordinate system.For example, consider a signal transduction mechanism 6-10 , hierarchical self-assembly 11-20 , Kinesin motor protein translocation on a microtubule 21, 22 , or hydrogen combustion H 2 /O 2 23, 24 .These systems transition from one stable configuration to another on the occurrence of some trigger event which may comprise of an external stimulus or the system reaching a special configuration.An external stimulus could be an input of energy that initiates hydrogen combustion, leading to a larger release of energy by the reaction itself.A special configuration could be a signaling molecule binding to an active receptor site.These stable configurations can be considered as fixed points in a transformed (reaction) coordinate system.The fixed points, the events, and their associated transitions are the coarse-level descriptions that are to be captured in the approximate N-body potential.However, it is still a challenge to construct a N-body Hamiltonian potential in a systematic manner that encodes the known coarse-level behaviors into a mathematical formulation and successfully predicts intermediate-scale transition events.This article introduces a method of encoding coarse-level dynamical behavior into logic functions that are used to "stitch" arXiv:1605.05348v1[physics.chem-ph]30 Apr 2016 together pairwise interaction potentials into an N-body potential.In this method, the practitioner uses experimentally observed coarse-level behavior to derive logic tables that capture various rules of interaction in the system.The qualitative logic tables are turned into a collection of quantitative logic functions associated with pairwise interaction potentials.The logic functions are then turned into smooth encoding functions via a replacement procedure which in turn are used to modify the pairwise potentials.The effect of an encoding function multiplying a pairwise potential is to smoothly turn the potential on or off when a precise set of conditions are met.The combination of the modified potentials gives an N-body potential that approximates the true potential governing the system.Using experimentally observed data and quantum calculations (red), we extract coarse-grain behavior (orange) namely: interactions rules and pairwise interaction potentials.This information is used to obtain an N-body potential (blue) for the system by employing the proposed formalism (green).
The method generates a potential that respects what is currently known about the system; it is not claimed that this method results in the unique potential governing the real system.The method does this by leveraging the existing experimental data and the coarse-level behavior that can be derived from it.If more experimental data becomes available, the same procedure can be used to generate a new potential that better models the system.This is equivalent to a refinement of the logic functions and ultimately a refinement of the generated potential.The resulting potential can have a much smaller dimension than the true potential and still accurately capture the relevant physics.
This article begins with a motivating example which is used as an impetus for our modeling framework.In the Methodology section, we define the major components of the framework -logic functions, permissible logical operations, and the translation to the associated encoding functions -and specify how they combine with the pairwise potentials to define the approximate potential.The procedure is depicted in Fig. 1.
The procedure is applied to three examples of increasing complexity to showcase the modeling framework.The first is a simple model of an inhibitor molecule mechanism.It shows how one would go from known coarse-level behavior to an approximate global potential that captures that behavior by explicitly constructing the logic tables, the associated logic functions, and the smooth encoding functions.The inhibitor molecule mechanism has more complicated logic than the motivating example and more effectively demonstrates the modeling procedure.
The second example shows how to model a simple, kinetically controlled, bond breaking chemical reaction using this framework.It shows that bond breaking events, and more generally chemical reactions requiring activation energy, can be naturally modeled in the framework.The general procedure for modeling a bond breaking event and how to account for the activation energy is shown.Furthermore, the derived potential is used with LAMMPS 25 to numerically simulate the chemical reaction.By changing the relative dissociation energies, the reaction can be biased in a particular direction.
As opposed to our method, many force fields have trouble capturing bond breaking events 26 .An exception is the ReaxFF potential 2 that was developed to model reactions of hydrocarbons.The derivation of ReaxFF is based on using interatomic distances to compute the bond order between two atoms and then using the bond order to obtain the bond energy.Corrections to the bond order are dependent on the valency and the deviation of the uncorrected bond order of an atom with its valency.Corrections to the bond energy, in the form of energy penalties (e.g. for over-/under-coordination) are added to get the system energy.This is contrasted with our method where bond weakening and breaking is due to the encoding function which is derived from coarse-level observed behavior.
The final example is a simple model of DNA transcription.It is shown that after the binding of RNA polymerase to the promoter region we can sequentially add the complementary base nucleotides to the DNA strand that is to be transcribed.DNA transcription is a complex process involving the interaction of many different molecules 27, 28 .This example shows that we can model such a complex process with a relatively low-dimensional potential that captures the observed mesoscale behavior.To the authors' knowledge there is no other other current potential accomplishing this task.

Motivating Example
There are a number of examples in biology where chemical reactions occurring within a cell are initiated by some signal or stimulus, followed by an ordered sequence of biochemical reactions.Often the term signal transduction is used to refer to such processes.One such example is the epidermal growth factor (EGF) signaling 9, 10 .Motivated by this example, we construct a hypothetical system to demonstrate how the proposed formulation can be used to construct a Hamiltonian potential for it.Each of the species in this system can be modeled as a rod having two sites of interaction at the end points, namely atoms {1, 2} on A, atoms {3, 4} on B, and atoms {5, 6} on C (Fig. 2).Let us write the force field energy for this system.In general, it is composed of the bonded energies formed from the stretch, bending, and torsional terms, the non-bonded van der Waals and electrostatic terms, and the coupling terms 26 .We can split the potential as Φ (i, j) + higher order terms where Φ (i, j) denotes a pair-wise interaction potential between two atoms i and j.These 2-atom terms encompass the stretching, torsional, van der Waals, and electrostatic terms and the higher order terms include the bending energies and all the k-atom (k ≥ 3) coupling terms.
For this system, the bonded energy terms are composed of the stretch and torsional energies between the atom-atom pairs (1, 2), (3,4) and (5,6), which we can group into a term U b .Assume that the only non-negligible non-bonded energy terms are the two van der Waals interactions between atoms 2 and 3 and atoms 2 and 5, and the coupling term between atoms 2, 3, and 5. Denoting these three terms by Φ (2,3) , Φ (2,5) , and Φ (2,3,5) , respectively, we get the force field energy of the system as The inclusion of the 3-atom potential Φ (2,3,5) is required in order to capture the transition of the pair (2, 3) being in a stable (bounded) configuration when atom 5 is not present to being in an unstable (free) configuration in the presence of the signaling atom 5.
While in general it may be hard to get the correct forms for the coupling term and other higher-order terms in the expansion, and thus the full potential, we know from the above observations that the effect of the potentials Φ (2,5) and Φ (2,3,5) is to basically to turn off Φ (2,3) when 5 is close to 2. Rewriting the potential as this means that term in parentheses is approximately 1 whenever atoms 2 and 5 are far and approximately 0 whenever atoms 2 and 5 are close.Instead of attempting to find the exactly functional forms of Φ (2,5) and Φ (2,3,5) , we approximate the potential as where S (2,3) is an encoding function that acts as a switch turning Φ (2,3) on and off.In this example, the encoding function is only function the distance between atoms 2 and 5.The encoding function takes values between 0 and 1, it is approximately 0 when atoms 2 and 5 are close, and approximately 1 when atoms 2 and 5 are far; thus it encodes the logic of the coarse-level observed behavior of the system.It is an approximation of the other terms: In the rest of the article, we make this approximation idea (eq.( 1)) precise and derive approximate N-body potentials from simple pairwise interactions that respect observed coarse-level behavior.We give a systematic procedure to construct the encoding functions which allows us to handle systems with more complex interaction rules.We will demonstrate the procedure with three examples.We also use molecular dynamics simulations using the derived potentials that show we can accurately capture the relevant physics.
There are a few items to keep in mind as motivation for the abstract concepts to follow.The basic building blocks for the N-body potential are pairwise interaction potentials (denoted by Φ (2,3) and Φ (2,5) for the above example).The explicit form of these potentials can be inferred from the experimental data or ab initio reasoning.We approximate the effect of the un-modeled potentials by modifying the relevant pairwise potentials with an encoding function.The encoding function only turns the corresponding potential on and off.The functional form of the potential does not change; it is only scaled between 0 and 1.The logic contained in the encoding functions is obtained from experimental observations or ab initio simulations and the logic only depends on pairwise distances between particles, except the pairwise distance in the argument of the encoding function's corresponding potential function; e.g. the logic corresponding to Φ (2,3) cannot depend on the distance between atoms 2 and 3.

Methodology
It is assumed that there are M interacting entities in a domain D, where D ⊆ R n , for n = 1, 2, or 3.Each entity is modeled by a finite number of particles with constraint forces between the particles; the totality of these particles over all the entities are labeled from 1 to N. This allows us to treat point particles as well as rigid and and non-rigid bodies.The configuration space is C = D N .An element, x x x ∈ C , takes the form x x x = (x x x 1 , . . ., x x x N ), where x x x j ∈ D describes the position of particle j in the domain D.
The dynamics of the system is driven by a potential gradient and external forces.Specifically, the functional form of the dynamics is where is the approximating potential.The notation in this equation needs to be explained.The set I ⊆ {1, . . ., N} × {1, . . ., N} is the set that defines which pairs of particles interact via some potential function.For example, if (2, 3) is in I , then there exists a pairwise potential between particles 2 and 3.The multiplicity function m : I → N determines the repetition of distinct elements of I ; the pair M = (I , m) is called a multiset and an element p p p ∈ I appears m(p p p) times in M .We make the convention that the pairs (i, j) and ( j, i) define the the same element of the multiset.Often m(p p p) will be 1 for every element of I , however, non-unit values become important when a pair p p p ∈ I interacts through multiple different potentials, each with its own encoding function.Such a non-unit multiplicity functions is useful when modeling, for example, bond-breaking chemical reactions.Two atoms interact through one potential, the bond potential, when they form a stable pair; when this bond is broken, another potential is required to model the electron-electron repulsion between the atoms.The term m i denotes the mass of particle i and ∇ x x x i is the gradient operator in the configuration space C with respect to the the position x x x i of particle i.The term F i ( x x x) collects the external forces on particle i.The external force not only contains the forces derived from any external fields, such as electric or magnetic fields, but can also be used to include stochastic effects or boundary constraints.The function Φ p p p, j : C → R denotes the j th potential that the pair p p p = (p 1 , p 2 ) ∈ I interacts through.The index j is an element of the set {1, . . ., m(p p p)}; i.e., it is the j th instance that pair p p p appears in the multiset.For x x x = (x x x 1 , . . ., x x x N ), it takes the form where for every i ∈ {1, . . ., N}, the coordinate map π i : C → D is linear and extracts the i th element, x x x i , of a configuration vector x x x = (x x x 1 , . . ., x x x N ); the norm • denotes the normal Euclidean norm; and φ p p p, j : R + → R is the specific 1D pairwise interaction potential through which the j th copy of pair p p p interacts.This potential could, for example, be a Lennard-Jones potential; it can also be different for different interaction pairs.The form given for the potential shows that it is only a function of the distance between x x x p 1 and x x x p 2 .The function S p p p, j : C → R in equation ( 2) is the encoding function, which will be defined shortly, corresponding to the potential Φ p p p, j ; it is the encoding function corresponding to the j th potential that atom pair p p p interacts through.It encodes the coarse-level interaction rules and it is a function of pairwise distances between particles, except for the particle pair p p p to which it corresponds.That is, for the interaction pair p p p = (p 1 , p 2 ) ∈ I , the encoding function S p p p, j is not a function of the distance x x x p 1 − x x x p 2 .Since the encoding functions and potentials are functions of relative distances only, equation ( 2) defines a Hamiltonian system 29 when we neglect the forces F i ( x x x).
A majority of the rest of the paper develops the encoding functions and their properties and shows how one would go from coarse-level interaction rules to encoding functions using a few examples.It is assumed that the coarse-level, interaction rules and the interaction potentials φ j : R + → R are known.These come from analyzing experimental data or ab initio simulations and are thus application specific and beyond the scope of this article.Ultimately, the encoding function S p p p, j will be a smoothed version, which is made precise later, of a logic function L p p p, j : C → {0, 1}.The logic function L p p p, j will be the constructed from a finite number of logical operations applied to elementary logic functions from a Boolean algebra.More precisely, the logic function will be an element of the Boolean sub-algebra generated by elementary logic functions.Thus to define the logic functions, it is required to know the specific definitions of the logical operators AND, OR, and NOT (symbolically denoted, ∧, ∨, ¬) and what Boolean functions are used as the elementary logic functions.
Before getting into the specific definitions, we would like to discuss what properties the logic functions should have and the motivation behind them.As hinted at above, the logic function L p p p, j encodes the coarse-level interaction rules.Its only effect is to turn the corresponding potential Φ p p p, j on and off depending on specific conditions met by the current configuration of the system.This means that L p p p, j has the effect of multiplying Φ p p p, j by either 0 or 1.As in the motivating example above, the logic function is only a function of pairwise distances like x x x i − x x x k .However, we choose to restrict L p p p, j so that it is not a function of x x x p 1 − x x x p 2 .Finally, it would be undesirable for L p p p, j to oscillate wildly when perturbing x x x.Since L p p p, j takes only the values 0 and 1, it will have discontinuities in general.The first property we wish L p p p, j to possess is that the set of points of discontinuity is "small" in the sense that its complement is open and dense in D. This means that for most configurations x x x, if x x x is perturbed slightly to y y y, then L p p p, j ( y y y) = L p p p, j ( x x x).The second property we wish the logic function to possess is similar but not equivalent to the first.Using the notation, r p p p = x x x p 1 − x x x p 2 , we can consider the logic function L p p p, j to be a function of the pairwise distances r q q q for q q q = p p p. If we fix all the r q q q 's except one, we wish for the resulting 1D function to make only finitely many transitions between 0 and 1.With these goals in mind, precise definitions of logic functions and elementary logic functions can be given.
A function b : C → {0, 1} is called a Boolean function on C and the set of such all such functions on C is denoted as B. It is easy to see that the functions that are identically 1 and 0 on C are Boolean functions.On B, define for all f , g ∈ B the two binary logical operations AND (∧) and OR (∨) and the unary logical operation NOT (¬) by It can be shown that (B, ∨, ∧, 1, 0) is a Boolean algebra with addition and multiplication given by ∨ and ∧, respectively.These three logical operations will be applied to specific elements of the set of all Boolean functions B on C to generate a Boolean sub-algebra.The logic functions L p p p, j will be elements of this sub-algebra.
The elementary logic functions generating the Boolean sub-algebra are formed from elements of the class of proximity functions.A proximity function P : R + → {0, 1} is a non-increasing function.What form do such functions have and why are they used as the elementary logic functions?Every proximity function has the form P R (r) = χ [0,R) (r), for some R satisfying 0 ≤ R ≤ ∞.The function χ [0,R) : R + → R is the indicator function for the semi-open interval [0, R) which takes the value 1 if the argument satisfies 0 ≤ r < R and 0 otherwise.Note that the functions that are identically 1 or 0 are proximity functions.The elementary logic functions are defined as compositions of a proximity function with the coordinate functions π i from above.Specifically, the elementary logic function q q q,R for particle pair q q q = (q 1 , q 2 ) ∈ I and parameter 0 ≤ R ≤ ∞ has the form This function is 1 when x x x q 1 and x x x q 2 are closer than distance R and 0 when not.A logic function L p p p, j is generated by applying finitely many of the logical operations ∧, ∨, and ¬ to the elementary logic functions.Recalling from above the desired properties for the logic functions, each L p p p, j corresponding to Φ p p p, j is only formed from finitely many logical operations on elementary logic functions q q q,R where q q q = p p p. Thus L p p p, j is a function of possibly all the pairwise distances r q q q = x x x q 1 − x x x q 2 except r p p p = x x x p 1 − x x x p 2 .
As we desired, each logic function L p p p, j is continuous almost everywhere in C .This follows since it is composed from pairwise products and sums of elementary logic function, which themselves are continuous almost everywhere.Since the logic function takes the value 0 or 1, this implies that around each point of continuity, there is a neighborhood on which L p p p, j is constant, as was desired.The second desired property follows from the constructing the logic function using a finite number of elementary logic functions.
Once the logic function is specified, it must be translated into a smooth encoding function.Ideally, this would be accomplished via a convolution in the (N-dimensional) configuration space with a smooth, nonnegative summability kernel (see 30 for a definition).Analytically, this is intractable, and computationally, this is very expensive.Instead, we individually smooth each of the 1D elementary logic functions q q q,R ( x x x) in the expression for L p p p, j ( x x x).This is done by replacing the proximity function of q q q,R with a smoothed version.Again, this could be done via the convolution (now 1-dimenional) of each proximity function with a smooth, 1D summability kernel or, alternatively, by the replacement of each indicator function with a specific functional form.In this article, we choose the latter approach and replace each proximity function χ [0,R) (r) with a function of the form For example, if the logic function has the expression then the corresponding encoding function would be for some choices of parameters α 1 , α 2 , n 1 , and n 2 .The parameters α and n control how well h α,n approximates a proximity function.In particular, h α,n (0) = 1 for all nonnegative α and positive n.For α > 0, lim r→∞ h α,n (r) = 0 and it is strictly monotonically decreasing.Note that for any n, h 0,n (r) = 1 for all r ≥ 0. On the other hand, for a fixed α > 0, the transition from 1 to 0 becomes sharper as n increases (Fig. 3b).To match a specific indicator function χ [0,R) , choose α = 1/R.With this choice of α, the function satisfies h 1 R ,n (R) = 1/2 for any n and as n → ∞,

Examples
To show the entire process, starting from coarse, interaction rules and recovering the encoding function, we apply the method to three examples in increasing order of complexity.The first example is a model for an inhibitor molecule system and is used to exhibit the core methodology of the modeling framework.This system can be considered as an extension of the signaling molecule example above (Fig. 2 use of the multiplicity function m from the framework.It is shown that the bond dissociation energy is accurately captured in this framework.Numerical simulations show that (i) the system exhibits the same coarse-level behavior that was used to derive the potential and (ii) that biased chemical reactions are easily handled.The final example is a simple model for DNA transcription and is the most complicated of the three.This example shows that the logic, and hence potential, of real systems can be captured in the modeling framework in a straight-forward manner.

Simple inhibitor molecule mechanism
This example can be thought of as a simple model for the action of an inhibitor molecule in a plane.Consider the three interacting molecules in Fig. 4. The configuration space for this example is C = (R 2 ) 6 , where x x x ∈ C is written as x x x = (x x x 1 , . . ., x x x 6 ).The set of interacting pairs is I = {(2, 3), (2, 5), (3, 6)}.For this example, the multiplicity function m : I → N is identically 1.Thus we have the pairwise potentials Φ (2,3) , Φ (2,5) , and Φ (3,6) .It is assumed that these potentials are formed using a Morse potential (see ( 12 This behavior is captured in the logic functions L (2,3) , L (2,5) and L (3,6) .The logic function L (2,3) is 0, i.e., the potential Φ (2,3) is turned off, when either atom 5 is close to atom 2 or when atom 6 is close to atom 3.This is different from the motivational example which only turned off the potential if 2 and 5 were close and the bonds between 2 and 5 or 3 and 6 never formed.Additionally, if AC has formed (atoms 2 and 5 close), then BC cannot form, i.e., L (3,6) = 0 and Φ (3,6) is off.Similarly, BC has formed (atoms 3 and 6 close), then AC cannot form.Table 1 captures this logic.
We need to specify what is meant by "close".We assume that "close" is in this case is determined by experiments to mean being within the distances R (2,5) and R (3,6) , respectively.Thus, atoms 2 and 5 are close when the elementary logic function (2,5),R (2,5) evaluates to 1 and not close when it evaluates to 0. Using the table, L 2,3 ( x x x) corresponding to the interaction potential 7/17 ) can be written as equation (7).
The other logic functions are To turn the logic functions into an encoding function, replace each of the proximity functions, χ (0,R p p p ] , in ( 7) -( 9) with their smooth versions, h 1/R p p p ,n p p p (eq. ( 6)).The encoding function S (2,3) where α p p p = 1/R p p p , and similarly for the other two.The approximate N-body potential for this system is The original configuration space was 12-dimensional.However, ( 11) is 8-dimensional since it only depends on four particles (four unique particles making pairs in I ).Thus we were able to reduce the dimension of the configuration space and still capture the relevant physics.
Here, we are only interested in demonstrating the methodology qualitatively so we make the approximation that derivative of each encoding function is 0 almost everywhere (this would be the case if the logic functions were used in place of the encoding functions in (11).With this approximation the force only consists of terms of the form S p p p ( x x x)∇Φ p p p ( x x x).One realization of the inhibitor molecule system (11) simulated in LAMMPS 25 with this approximation is shown in Fig. 5.The potentials Φ p p p are formed from Morse potentials where D is the dissociation energy, r eq is the equilibrium distance of the bond, and a is a parameter.Simulations are performed by solving the Langevin equations at constant temperature (i.e.NVE ensemble).The parameters used in the computations are given in Supplementary Table I.Initially, the AB complex is formed.Around 30 femtoseconds C comes close enough, turns off the AB bond and BC forms and can diffuse away from A. This remains the case until around 450 fs, when A approaches BC, the BC bond turns off and the AC bond turns on.Supplementary Movie 1 shows a simulation of the inhibitor molecule mechanism.

Modeling a bond breaking chemical reaction
Let us consider a slightly more complicated example than the inhibitor molecule.We model a bond breaking, reversible, chemical reaction.In particular, we will model the reaction Modeling such chemical reactions is difficult with traditional force field methods since they cannot describe changes in the electronic structure and, thus, are unable to describe bond-breaking, bond-forming, charge transfer, etc., of the system undergoing a reaction 4, 31 .Rather than solving the quantum mechanical equations, we take a coarse-level approach and approximate the bond breaking mechanism with logic functions.This example makes use of the multiplicity function m(i i i) in ( 2) in order to model the electron-electron repulsion during the transition state.It also shows that the use of a smooth encoding function accurately accounts for the bond dissociation energy.Let x x x = (x x x 1 , . . ., x x x 6 ) ∈ (R 2 ) 6 be the configuration vector for this system (see Fig. 6).Table 3. Bond-breaking logic rules.
The logic functions are logical NOT's of x x x 2 -x x x 5 and x x x 1 -x x x 3 proximity functions: We assume that r eq AB < R dis AC,B and r eq AC < R dis AB,C .The use of the smooth encoding function in the potential (as opposed to the logic function) allows the transfer of the correct amount of energy from C to AB in order to break the bond; C must transfer an amount of energy equivalent to the bond dissociation energy D AB of the AB bond in order to turn off the Φ (1,3),1 potential.We refer the reader to Sec.III.A in the Supplementary Information for the derivation.
Let us now turn to the case occurring directly after a successful collision of C with AB.In this case, A and B are close to their equilibrium distance ( x x x 1 − x x x 3 ≈ r eq AB ) and A and C are closer than the AB-bond dissociation distance ( x x x 2 − x x x 5 < R dis AB,C ).In this state, the bonds are weak and neither AB nor AC is stable; the system is at its transition state.In this transition state the forces experienced by the molecules due to the bond potentials Φ (2,3),1 and Φ (2,5),1 are small since the encoding functions and their partial derivatives are small, and thus the bond potentials are approximately "off".The dynamics are predominantly dominated by noise and the residual momentum of the molecules.
In this transition state, the electron-electron repulsion should be directly accounted for via a short-range repulsion potential between the molecules, such as by using the repulsive part of a Morse potential.The logic functions are defined such that these repulsion forces are only "on" when the system is in its transition state.This is easily accomplished.Denote the short-range repulsion potential between A and C by Φ (2,5),2 .This force is defined such that Φ (2,5),2 ( x x x) ≈ 0, for x x x 2 − x x x 5 > R dis AB,C .This force is turned on when x x x 1 − x x x 3 < R dis AC,B .The logic function for the A-C repulsion is Similarly, repulsions between A-B and B-C can be defined with logic functions similar to the above one.There are three possible outcomes for when the system exits its transition state: (i) either AC forms a stable molecule, (ii) AB reforms, or (iii) no bonds are formed and all the molecules are free molecules.This depends on the equilibrium distances of the bonds, the dissociation distances, the incoming momentum of C, and the repulsion forces.Figure 7 shows the two most probable outcomes for a single AB + C event.
For simulations, the AB and AC bonds (Φ (1,3),1 and Φ (2,5),1 , respectively) are given by Morse potentials, (12).In simulations, only the short-range A-B and A-C electron-electron repulsions are modeled and are only active during the transition state.The form of these for these repulsions are chosen as the repulsive part of a Morse potential with the same parameters as the full potentials used for the AB and AC bonds.The logic function for the A-C repulsion potential, Φ (2,5),2 , is given by ( 14) with obvious modifications for Φ (2,5),2 .The associated encoding functions are given by the normal replacement procedure.The full potential used during the numerical experiments is given in (15).
D AC e −2a( x x x 2 −x x x 5 −r eq AB ) − 2e −a( x x x 2 −x x x 5 −r eq AC ) D AB e −2a( x x x 1 −x x x 3 −r eq AB ) D AC e −2a( x x x 2 −x x x 5 −r eq AB ) The force derived from (15) is used in LAMMPS 25 to simulate the system for an unbiased and a biased potential (parameters in Supplementary Information Table II).The parameters of the first simulation are chosen so that the AB and AC are symmetric (D AC /D AB = 1).In this case, the chemical reaction is unbiased and if averaged over all realizations of the noise, it is expected that the amount of time AB is formed is equal to the amount of time AC is formed.Figure 8 shows the potential energy for this simulation (Fig. 8a), the corresponding level sets (Fig. 8b), and a typical realization of the simulation (Fig. 8c).In the energy surface plot and the level set plot, the symmetry of the potential is evident.The realization shown in Fig. 8c starts with AB near its equilibrium length (2 Å) with C far from A. The realization shows the approximately equal times that AB and AC are formed.The deviation is due to this being a particular realization rather than an average over an ensemble of realizations and the finite nature of the simulation.
The parameters of the second simulation are chosen so that the reaction is biased in favor of AC.With the chosen parameters (D AC /D AB = 2), the AC is twice as stable as AB.If the noise in the system is strong enough, it is expected that AC is formed twice as often as AB is formed.Figure 9 shows the potential energy for this simulation (Fig. 9a), its corresponding level sets (Fig. 9b), and a realization of the simulation (Fig. 9c).In the energy surface plot and the level set plot, the symmetry of the potential is evident.The realization shown in Fig. 9c starts with AB near its equilibrium length (2 Å) with C far from A. In this particular realization AC forms very quickly.Figure 9c shows the bias towards the more stable AC.The system spends most of its time with a stable AC molecule with a relatively small amount of time with a stable AB molecule.Thus, biased reactions can be captured in the framework.A movie of a part of the unbiased reaction simulation can be found in Supplementary Movie 2.

DNA transcription model
The final example is inspired by DNA transcription 27, 28 .The model consists of a promoter region (sites 1 & 2) to which RNA polymerase (RNA pol) binds (sites 3 & 4), and a four nucleotide DNA strand, ACTG, to be transcribed (Fig. 10).
12/17  As a first approximation of the transcription process, the movement of the RNA polymerase down the DNA chain and the unwinding/rewinding of the DNA have not been explicitly modeled.
In the absence of the RNA pol, the free nucleotides cannot bind to their complementary nucleotides in the 4 nucleotide DNA strand (ACTG = (5,6,7,8)).Once RNA pol binds to the promoter, the first first nucleotide (A) in the DNA strand can bind to the free version of its complementary nucleotide (U).Before this binding happens, the remaining nucleotides in the strand (CTG) cannot bind with their (free) complementary nucleotides.Once A has bound to a free U nucleotide, the next nucleotide in the strand (C) can bind with a free G nucleotide, while the remaining two nucleotides (TG) still cannot bind with their complementary nucleotides.Once the free G has bound with C, the sugar and phosphate groups on T and G can bind to start forming the backbone of the complementary strand.This cascading process continues until each nucleotide in the original DNA strand ACTG has bound with its complementary nucleotide, resulting in the complementary RNA strand UGAC.At this point, the complementary strand and the RNA pol unbind from the original strand and promoter region, respectively.
Supplementary Table III lists the reaction potentials for each of the interacting pairs.The nucleotide base pairs interact via a hydrogen bond φ H , whereas the sugar and phosphate groups covalently bond through φ SP .The interaction potentials for the system can be easily read from this table (see Supplementary Information Sec.IV.A).
Let us step through the rest of the logic in the order the reaction occurs: (i) The bonds between the RNA pol and the promoter region (Φ (1,3) and Φ (2,4) ) are "on" except when the complementary chain has formed and is still attached to the original base strand.(ii) The A-U bond (Φ (5,9) ) is "on" when the RNA pol has bonded with the promoter and the complementary chain has not formed.This second condition prevents the complementary strand from reattaching to the original DNA strand once it has been formed.It is "off" otherwise.(iii) The C-G bond (Φ (6,12) ) is "on" when all of following conditions are true: (1) RNA pol has bonded with the promotor, (2) the A-U bond has formed.It is "off" otherwise.(iv) The sugar-phosphate group bond Φ (11,13) turns "on" after the C-G bond has formed.It is "off" otherwise.(v) The T-A bond (Φ (7,15) ) is "on" when all of the following conditions are true: (1) RNA pol has bonded with the promotor, (2) the A-U bond has formed, (3) the C-G bond has formed, and ( 4) the (11, 13) sugar-phosphate bond is formed.It is "off" otherwise.(vi) The sugar-phosphate group bond Φ (14,16) turns "on" when T-A bond has formed.It is "off" otherwise.(vii) The G-C bond (Φ (8,18) ) is "on" when all of the following conditions are true: the conditions in (v) are true, the T-A bond has formed, and the (14,16) sugar phosphate bond has formed.It is "off" otherwise.(viii) The sugar-phosphate group bond Φ (17,19) turns "on" when G-C bond has formed.It is "off" otherwise.
A global potential derived from the above logic rules is given in Eq. ( 16).The exact form of the logic functions L p p p ( x x x) and the associated smooth encoding functions S p p p ( x x x) comprising the potential are given in the Supplementary Information.The derivation of the potential is not difficult, but lengthy.We refer the reader the Supplementary Information Sec. 4 for the details.
sugar-phosphate backbone for complementary RNA strand (16) Figure 11 shows a trace of the pairwise distances between atoms for a typical simulation using this potential in LAMMPS (parameters in Supplementary Information Table IV).We use the same qualitative approximation of the force as was used in the inhibitor molecule example.For simplicity, all the potentials are taken to be Morse potentials.The red trace (TF) corresponds to the distance between the RNA pol and the promoter region.The variables r 5,9 (teal), r 6,12 (black), r 7,15 (orange), and r 8,18 (blue) correspond to the sites on the complementary A-U, C-G, T-A, and G-C pairs from the base strand and the free nucleotides.At the start, the RNA pol and the free nucleotides diffuse around in space.Around 900 ps, the RNA pol binds to the promoter region (TF trace ≈ 0).The free nucleotides then bind in the the order of the designed logic.U binds to A (r 5,9 ≈ 0) around 1100 ps; G binds with C (r 6,12 ≈ 0) between 1800 and 1900 ps; A binds to T (r 7,15 ≈ 0) around 2400 ps; and finally C binds to G (r 8,18 ≈ 0) around 2700 ps.Once this final free nucleotide has bounded with its complement, the complementary chain has finished forming and unbinds as does the RNA pol.The RNA pol can rebind to the promoter region, but the complementary RNA strand cannot rebind to the original DNA strand.This is exactly the behavior designed into the potential.Supplementary Movie 3 in the Supplementary Information shows one simulation of the DNA transcription.Around 900 ps, the RNA pol binds to the promoter region (TF trace ≈ 0).The free nucleotides then bind in the the order of the designed logic.U binds to A (r 5,9 ≈ 0) around 1100 ps; G binds with C (r 6,12 ≈ 0) between 1800 and 1900 ps; A binds T (r 7,15 ≈ 0) around 2400 ps; and finally C binds to G (r 8,18 ≈ 0) around 2700 ps.Once this final free nucleotide has bounded with its complement, the complementary chain is finished formed and unbinds as does the RNA pol.The RNA pol can rebind to the promoter region, but the complementary strand cannot rebind to the original DNA strand.

Conclusions
We have developed and demonstrated a methodology and mathematical framework for obtaining an approximate interaction potential for a system which respects known coarse-level behavior.This methodology develops a semi-empirical model for the system by encoding the known coarse-level physics into logic functions that then modify simple pairwise potentials.Each logic function's only role is to turn its associated pairwise potential on or off.A smooth multi-body interaction potential is obtained by replacing each logic function with a smoothed variant.The reader may wish to think of the resulting approximate potential as a linear combination of pairwise potentials where instead of the coefficients taking scalar values, they are encoding functions capturing the coarse-level logic.Three relatively simple examples demonstrated our methodology: a simple inhibitor molecule mechanism, a chemical reaction with bond breaking, and a model inspired by DNA transcription.While these examples were simple and inspired by biophysical and chemistry applications, we stress that the methodology is quite general and not restricted to these application domains or only simple problems.Any system that is driven by a potential can utilize this methodology to its benefit.
The result of our procedure is the approximation of a complicated, high-dimensional potential with a lower-dimensional representation that still respects the relevant physics.A significant reduction in the dimensionality of the system is possible; instead of accounting for every interaction between a large number of components, we now only need as many variables as are needed to correctly model the coarse-level logic.In the bond breaking example, the potential capturing the logic was 8-dimensional, whereas the dimension of the configurations space was 12.The same system modeled at the quantum level is much more complicated.Since the bond breaking event is the relevant physics, the reduced order model is accurate enough for this purpose.
With this dimensional reduction, the ability to accurately simulate large, complicated systems within a computational design framework is feasible.The resultant models can be wrapped in an optimization loop as part of exploratory computational experiments, such as for the development of new drug therapies, or as part of an engineering design loop.This in turn allows for the faster and cheaper development of new technologies and products.
We note that the developed framework can be potentially used in reverse: not for approximation to a given physical process with coarse-grained logic given, but for design of molecular processes with logic prescribed by a designer.This is achieved by providing to the designer the specifications of molecules that can carry the logic out.

Figure 1 .
Figure 1.Flow chart of procedure.Using experimentally observed data and quantum calculations (red), we extract coarse-grain behavior (orange) namely: interactions rules and pairwise interaction potentials.This information is used to obtain an N-body potential (blue) for the system by employing the proposed formalism (green).

Figure 2 .
Figure 2. Simple signaling molecule mechanism, A + B → AB C − → A + B, modeled by three rods.When C is not present, A and B form a complex.When C is present, A and B dissociate and diffuse apart.Molecule C is free to diffuse away from molecule A. This behavior is captured with the following rules.When atom 5 is far from atom 2, the potential between atoms 2 and 3 is on.When atom 5 is close to atom 2, the potential between atoms 2 and 3 is turned off allowing molecules A and B dissociate and diffuse apart.Atom 5 can diffuse away from atom 2.

Figure 3 .
Figure 3. Behavior of the h α,n function from equation (6)

Figure 4 .
Figure 4. Inhibitor molecule example.When the inhibitor molecule, C, is not present, a bond between receptor A's site 2 and site 3 on the active molecule B can form.When the inhibitor molecule is present, it can either take up the receptor site through a (2,5) bond or bind to site 3 on B with site 6.Either of these cases prevents to active molecule B from binding with its receptor site on A.

Figure 5 .
Figure 5.One realization of the inhibitor molecule system (11) simulated in LAMMPS.Initially, the AB complex is formed.Around 30 femtoseconds C comes close enough, turns off the AB bond and BC forms and can diffuse away from A. This remains the case until around 450 femtosceonds, when A approaches BC, the BC bond turns off and the AC bond turns on.

Figure 6 .
Figure 6.Diagram for chemical reaction example.Φ (1,3),1 and Φ (2,5),1 represent the stable bonds AB and AC, respectively.The dashed lines represent the repulsion forces induced by the encoding functions.In particular, in the left diagram, the repulsion between A and C is due to the partial derivatives of S (1,3),1 ( x x x) with respect to both x x x 2 and x x x 5 .See Supplementary Information Sec.III.A for a discussion of the repulsion force induced by the smooth encoding functions.

Figure 7 .
Figure 7.The two most probable outcomes of a successful AB + C event.Both trajectories start at the same configuration the difference is that C has a greater momentum for the blue (thicker) curved arrow.For both the red and blue trajectories, C approaches A. The AB bond breaks when r AC = x x x 2 − x x x 5 < R dis AB,C .Depending on the momentum and the relative strength of the repulsion terms, either AB reforms or AC forms.

( a )Figure 8 .
Figure 8. Simulation of an unbiased (1:1 well-depth), bond breaking chemical reaction, (13).(a) The potential energy (15) for the system.The parameters are given under simulation 1 in Supplementary Table II.(b) The level set plot of the potential energy.(c) A typical trajectory of the simulation.The cyan trace denotes the distance between molecules A and C (r AC = x x x 2 − x x x 5 ), whereas the red trace corresponds to the distance between molecules A and B (r AB = x x x 1 − x x x 3 ).Initially, A and C are near their equilibrium length (2 Å) and B is far from A. We see a successful AC + B → AB + C event happening very soon (red trace is close to the equilibrium distance, then becomes large; cyan trace is large, then becomes small).

( a )Figure 9 .
Figure 9. Simulation of biased (2:1 well-depth), bond breaking chemical reaction, (13).(a) The potential energy (15) for the system.The parameters are given under simulation 2 in Supplementary Table II.(b) The level set plot of the potential energy.(c)A typical trajectory of the simulation.The cyan trace denotes the distance between molecules A and C (r AC = x x x 2 − x x x 5 ), whereas the red trace corresponds to the distance between molecules A and B (r AB = x x x 1 − x x x 3 ).Initially, A and B are near their equilibrium distance (2 Å) and C is far from A. We see a successful AB + C → AC + B event happening very soon (red trace is close to the equilibrium distance, then becomes large; cyan trace is large, then becomes small).The trace exhibits the bias towards a stable AC bond, since the cyan trace is close to equilibrium longer than the red trace.

Figure 10 .
Figure 10.Simple DNA transcription model.Free base nucleotides cannot bind with the DNA strand until RNA pol binds with the promoter.When RNA pol is bound to the promoter, the free nucleotides bind to the DNA strand ACTG sequentially from left to right.Once the complementary strand has formed, the RNA pol unbinds from the promoter and then the complementary strand can diffuse away.Once the RNA pol has diffused far enough away, the bonds between the complementary base pairs turn off and the complementary strand can diffuse away.Dashed arrows between sites denotes an active potential.The P blocks denote a phosphate group and the S blocks denote a sugar group.

Figure 11 .
Figure11.DNA transcription.The red trace (TF) corresponds to the distance between the RNA pol and the promoter region.r 5,9 (teal), r 6,12 (black), r 7,15 (orange), and r 8,18 (blue) correspond to the sites on the complementary A-U, C-G, T-A, and G-C pairs from the base strand and the free nucleotides.At the start, the RNA pol and the free nucleotides diffuse around in space.Around 900 ps, the RNA pol binds to the promoter region (TF trace ≈ 0).The free nucleotides then bind in the the order of the designed logic.U binds to A (r 5,9 ≈ 0) around 1100 ps; G binds with C (r 6,12 ≈ 0) between 1800 and 1900 ps; A binds T (r 7,15 ≈ 0) around 2400 ps; and finally C binds to G (r 8,18 ≈ 0) around 2700 ps.Once this final free nucleotide has bounded with its complement, the complementary chain is finished formed and unbinds as does the RNA pol.The RNA pol can rebind to the promoter region, but the complementary strand cannot rebind to the original DNA strand.

Table 2 .
Potentials in chemical reaction model.