Many biomolecular systems or other complex macromolecules can dynamically visit a broad range of conformational states. External perturbations such as a molecular interaction or a mechanical force can cause a molecule to dynamically transit between these conformational states. While the conformational space of biomolecules is typically analyzed by coordinate-based methods such as the detection of correlated motions, Force Distribution Analysis (FDA) has been recently developed as an alternative approach to analyze structure and structural transitions [1]. The advantages of analyzing internal forces instead of coordinates are two-fold. First, forces between atoms or residues represent internal coordinates, which consequently do not require any structural fitting. Secondly, forces are a more sensitive measure, as they are able to reveal low-amplitude yet functionally important motions such as those in a stiff protein core [2]. Among others, internal forces obtained from FDA proved able to explain the mechanical robustness of immunoglobulin domains [3] and to reveal the pre-stress in ubiquitin [4].

FDA has been applied so far to averaged dynamical data from Molecular Dynamics (MD) simulations. However, the dynamics of the force distribution within proteins or other macromolecules, *e.g.* in equilibrium, under a constant load, a load varying in time, or upon binding of another molecule, can only be characterized by following the changes of the internal forces in time, which cannot be easily achieved with the previous implementation of FDA. We here introduce a Time-Resolved Force Distribution Analysis (TRFDA) method, which adds a temporal component to FDA to enable the analysis of pairwise forces associated with conformational changes.

### Implementation

#### Overview

Time-Resolved Force Distribution Analysis (TRFDA) is based on the same concept of using pairwise forces as FDA, but focuses on their evolution during the MD simulation as well as on the analysis of large molecular systems. To achieve these goals, it stores in memory pairwise forces from only one integration time step at a time, such that the memory usage is independent from the length of the simulation. It also saves output data in new file formats, which are well suited for time series analysis. Similar to FDA, atomic pairwise forces are computed for all types of bonded interactions as well as from Coulomb and Lennard-Jones potentials; long-range electrostatic interactions computed on a grid (*e.g.* Particle Mesh Ewald [5, 6]) cannot be decomposed in atomic pairwise forces. As in the previous FDA implementation, an MD simulation is first performed to obtain a trajectory, which is then used as the basis for a rerun using the TRFDA code. The core code is completely rewritten for GROMACS 4.5.3 [7], and VMD [8] plugins are provided for visualization. The implementation includes several notable improvements over FDA, which are discussed below.

#### Decomposition of forces from 3- and 4-body potentials

While two-body potentials such as bonds, Coulomb and Lennard-Jones potentials, can be directly used for the analysis of pairwise forces, many-body potentials which act on more than two atoms need to be decomposed into pairwise forces. TRFDA introduces a complete decomposition of the forces resulting from 3- and 4-body potentials (angle, dihedral angle, cross bond-bond, cross bond-angle), as described in the Additional file 1. The pairwise forces resulting from this way of force decomposition can have any direction, *i.e.* do not align with the vector connecting the two atoms involved. To accommodate this, all pairwise forces are stored and handled internally as vectors and are transformed into scalar values only when written to a file, if desired by the user.

The decomposition rectifies the shortcomings of the previous FDA code, which used approximations for computing pairwise forces resulting from 3- or 4-body potentials. For an angle formed by atoms *i*, *j* and *k* (*j* between *i* and *k*), it considered no pairwise forces to act on atom *j*, and assumed the pairwise force between atoms *i* and *k* to act along the distance vector between the two atoms. A similar decomposition has been used for a dihedral angle formed by four atoms. In contrast to these approximations, the decomposition used in TRFDA correctly reproduces the distribution of angle and dihedral forces in a molecule.

#### Internally computed pairwise forces between residues

Computing pairwise forces between residues allows a significant decrease in storage and computational cost for analysis, while providing a mapping of the interactions on the primary and possibly also secondary structure of a protein. The examination of pairwise forces between residues can also be used as a tool in the development of residue-level coarse grained models [

9]. The pairwise force representing the interaction between residues

*ri* and

*rj* acts on the centers of mass of the 2 residues and is calculated as:

where *i* is an atom of residue *ri* and *j* is an atom of residue *rj*, with *ri* and *rj* being different. TRFDA computes internally pairwise forces between residues. As a consequence, TRFDA has relatively small memory requirements and completely avoids writing out pairwise forces between atoms, if only forces between residues are of interest. This is a significant improvement over the previous implementation of FDA, which involved storing in memory a large number of pairwise forces between atoms and saving them in a large file, which is subsequently read by a standalone tool.

A further advantage of TRFDA is the equal treatment of atoms and residues with respect to the output options. For example, the same vector to scalar transformations can be applied to both pairwise forces between atoms and pairwise forces between residues, punctual stress (see below) can be calculated both per atom and per residue, and the output file formats are the same for data referring to atoms or to residues.

#### Summed versus detailed pairwise forces

In TRFDA, the force decomposition described above allows several pairwise forces to be computed for the same atom pair. For example, two atoms forming a bond might also form angles and dihedral angles with neighboring atoms, resulting in multiple forces between them from different potentials. For two atoms, it is most often interesting to calculate the total interaction between them, independent of the underlying energy potentials. Only few applications, like force field development, might require a differentiation based on the potentials, with an associated increase in memory usage and output file size. TRFDA accommodates both scenarios: in the first case, it computes a vector sum from all pairwise forces between the two atoms, while in the second case it computes separate vector sums for pairwise forces resulting from different types of potentials.

In contrast, due to the approximations used for 3- and 4-body interactions in the previous FDA version, the bonded pairwise force between two atoms results from exactly one type of interaction.

#### Choice of output

TRFDA comes with a wide range of output options. Although internally forces are calculated and stored only in vector form, TRFDA can output pairwise forces as vector or scalar values, and several derived quantities based on them.

A scalar pairwise force is computed as the magnitude of the vector pairwise force, same as in FDA, or as the magnitude of the pairwise force projected onto the distance vector between the two atoms, same as in [4]. In either case, the scalar value carries a sign indicating whether it is a repulsive (plus) or attractive (minus) force. If the angle between the vector pairwise force and the distance vector between the atoms is in the range (- *Π*/2, *Π*/2), the pairwise force is considered attractive. If the angle between the vector pairwise force and the distance vector between the atoms is in the range (*Π*/2, 3 *Π*/2), the pairwise force is considered repulsive. If the vector pairwise force is perpendicular to the distance vector, the pairwise force can be considered neither attractive nor repulsive and it is set to zero.

From the absolute values of the scalar pairwise forces acting on an atom, TRFDA can compute the sum and average as well as select the minimum or maximum. The sum of the absolute values of scalar pairwise forces acting on an atom

*i*:

measures the stress acting on that atom. It is reminiscent of the calculation of stress from continuum models used in mechanical engineering. We denote the stress defined in Eq. 2 as *punctual stress*, emphasizing the action of the force on a dimensionless point instead of an area, which is here ill-defined. As a sum of forces, the punctual stress is expressed in units of force. Such a deviation from the typical definition of stress is also found for the atomic virial stress [10], expressed in units of energy, and stems from the difficulty of defining geometrical properties (area for punctual stress or volume for virial stress) at atomic level. For systems in equilibrium, the sum of scalar pairwise forces from Eq. 2 can converge to non-zero averages over time; in contrast, the sum of vector pairwise forces acting on a single atom *i*, typically computed in MD simulations, averages to zero in the same conditions. As opposed to the analysis of pairwise forces, which gives a very detailed view of their distribution in or between molecules, the punctual stress expresses in a simple way where pairwise forces accumulate and allows the detection of atomic-level “hot-spots”.

For both pairwise forces and per atom quantities, the output consists of simple to parse text files. Their format allows an unlimited number of atoms, pairwise forces and simulation time steps, and is documented in the manual, which accompanies the code. We chose text as opposed to binary file formats to allow for loading the data in various analysis tools at the expense of the size of the files.

TRFDA can also write out scalar pairwise forces in the same format as FDA, preserving compatibility with the FDA R library [1]; however the results should be interpreted with care, as the 3- and 4-body interactions are expressed differently.

#### Internal organization

TRFDA defines two groups of atoms between which pairwise forces are computed. A pairwise force is only computed when the two atoms belong to different groups, while pairwise forces between atoms of the same group are not computed. This allows, for example, interactions between a protein and a ligand to be computed efficiently. If the interactions inside the protein and ligand are also of interest, the two groups should comprise the same atoms, namely those of both protein and ligand. This way of selection avoids the calculation of possibly unnecessary pairwise forces within parts of the molecules, in contrast to the previous implementation of FDA [1].

For internal storage of pairwise forces, FDA uses a square matrix-like structure for which the memory usage depends on the atom numbering in the molecular system. TRFDA replaces this structure with lists for a much more efficient memory usage, opening the possibility of calculating pairwise forces for significantly larger molecular systems. Using lists also allows storing, for the same two atoms, of separate pairwise forces from the different potentials, *i.e.* the detailed pairwise forces described above. Searching in lists is not as fast as the direct access in the matrix-like structure, but this is partly offset by the better locality of data access.