11institutetext: French Council of Scientific Research (CNRS), LIRIS, Lyon,
France
11email<EMAIL_ADDRESS>
22institutetext: University of Lyon, ERIC-LAFMIA, Bron, France
22email<EMAIL_ADDRESS>
# Building Analytics Pipelines for Querying Big Streams and Data Histories
with H-STREAM ††thanks: The authors would like to thank Enrique J. Arriaga-
Varela who participated in the implementation of the first version of H-STREAM
funded by the Mexican CONACyT Young Scientists’ Support Program. Experiments
were done at the Extreme Computing group of the Barcelona Supercomputing
Centre. We thank Microsoft for the Azure grant we obtained that enabled the
scalability tests we performed.

Extended Version

Genoveva Vargas-Solar 11 Javier A. Espinosa-Oviedo 22

###### Abstract

This paper introduces H-STREAM, a big stream/data processing pipelines
evaluation engine that proposes stream processing operators as micro-services
to support the analysis and visualisation of Big Data streams stemming from
IoT (Internet of Things) environments. H-STREAM micro-services combine stream
processing and data storage techniques tuned depending on the number of things
producing streams, the pace at which they produce them, and the physical
computing resources available for processing them online and delivering them
to consumers. H-STREAM delivers stream processing and visualisation micro-
services installed in a cloud environment. Micro-services can be composed for
implementing specific stream aggregation analysis pipelines as queries. The
paper presents an experimental validation using Microsoft Azure as a
deployment environment for testing the capacity of H-STREAM for dealing with
velocity and volume challenges in an (i) a neuroscience experiment and (in) a
social connectivity analysis scenario running on IoT farms.

###### Keywords:

Stream processing Cloud Micro-Services.

##  1 Introduction

The Internet of Things (IoT) is the network of physical devices, vehicles,
home appliances, and other items embedded with electronics, software, sensors,
actuators, and connectivity, enabling these objects to connect and exchange
data. When IoT is augmented with sensors and actuators, it enables the
construction of cyber-physical systems, including smart grids, virtual power
plants, smart homes, intelligent transportation and smart cities. These
systems are data centred where data are produced at different paces and
requires an agile and resources aware management to be consumed, archived and
analysed, optimising in particular underlying computing, storage and
communication resources. Streams processing calls for research that can
intelligently process raw data and perform ”online” analysis combining
streaming data and historical persistent data, that could been collected by
IoT devices and that has been stored for postmortem analysis. Existing
platforms provide efficient solutions for processing streams with parallel
execution backends for example Apache Flink 111https://flink.apache.org, Kafka
222https://kafka.apache.org/intro and message-based infrastructures like
Rabbit MQ 333https://www.rabbitmq.com. Programmers rely on these platforms to
define stream processing operations that consume ”mini”-batches of streams
observed through temporal windows. The processing operations can be rather
complex and ad-hoc to target application requirements. It is technically
challenging for analytics-based applications to analyse data postmortem
combining it with streams to perform more agile and online analytics tasks.
For example, compare the number of people entering a shopping mall in the last
hour with the average number of people entering the shopping mall at the time
interval every day during the last month. This task is, of course, possible to
compute but requires data management and programming skills for doing it
properly.

Current advances in data processing and data analytics have shown that it is
possible to propose general operations as functions or operators that can be
called, similar to queries within databases applications (e.g. Spark
programs). Still, the stream/data processing operations remain embedded within
programs, and this approach prevents reutilisation and can imply high
maintenance costs. Providing platforms that can let combine processing, and
analytics operators and apply them on streams and stored data, in the spirit
of AI Microsoft Gallery are still to come. The challenge is to define an
architecture/platform with the right components that can wrap operators,
manage data collections and streams and execute pipelines of analytics tasks
that can continuously deliver analytics results to be consumed by target
applications.

This paper proposes H-STREAM 444https://youtu.be/WmTgGHY_FUY a big stream/data
processing pipeline enactment engine. It provides stream processing micro-
services555Micro-service is a software development technique that structures
an application as a collection of loosely coupled services. In a micro-service
oriented architecture, services are fine-grained and the protocols are
lightweight. for supporting the analysis and exploration of big streams in IoT
environments. H-STREAM combines stream processing and data storage techniques
that are tuned depending on the number of things producing streams, the pace
at which they produce them, and the physical computing resources available for
processing them on-line and delivering them to consumers. H-STREAM deploys big
stream operators pipelines, called big stream queries, on message queues
(Rabbit MQ) and big data processing platforms (Spark) to provide a powerful
execution environment.

The experiments present stream processing micro-services composed and
installed in a cloud environment for implementing specific stream aggregation
queries. The paper particularly focuses on two experiments using Microsoft
Azure as deployment environment introducing respectively velocity and volume
challenges of an IoT farm.

Accordingly, the remainder of the paper is organised as follows. Section 2
introduces related work regarding stream processing and data harvesting
looking into IoT environments. The section discusses some limitations and
underlines how our work intends to go a step forward concerning those
solutions. Section 3 describes the general architecture of H-STREAM with
operators as micro-services that are deployed on high performance underlying
infrastructures. It also introduces the micro-services composition language it
proposes for defining Big Stream processing pipelines as queries. Section 4
introduces the core of our contribution, a stream processing micro-service; it
describes its architecture and the data processing workflow it defines.
Section 5 describes the experimental scenario that applies series of micro-
services for processing streams produced continuously for analysing connection
logs to evaluate scale up in terms of the number of things and volume of
streams. Section 6 concludes the paper.

##  2 Related work

We analyse stream management systems that emerged to query streams from
continuous data providers (e.g. sensors, things). Stream processing is the
processing of data in motion, or in other words, computing on data directly as
it is produced or received. Big Streams can result from a fine-grained
continuous reading of the environment, society, and organisations, natural and
social phenomena. Observations are done in different conditions and with
different devices, and therefore, Big Streams are raw heterogeneous
continuously produced data flows that must be processed for extracting useful
information. The systems that receive and send the data streams and execute
the application or analytics logic are called stream processors. The primary
responsibilities of a stream processor are to ensure that data flows
efficiently and the computation scales and is fault-tolerant.

Many stream processors adopt stateful strategies. Stateful stream processing
[5, 4, 2] maintains contextual state used to store information derived from
the previously-seen events. Before stream processing, streams were often
stored in a database, a file system, or other forms of mass storage.
Applications would query the data or compute over the data as needed.

Apache Storm666https://storm.apache.org is a distributed stream processing
computation framework that is distributed, fault-tolerant and guarantees data
processing. A Storm application is designed as a ”topology” in the shape of a
directed acyclic graph (DAG) with spouts and bolts acting as the graph
vertices. Edges on the graph are named streams and direct data from one node
to another. Together, the topology acts as a data transformation pipeline.
Apache Flink777https://flink.apache.org/ is an open-source stateful stream
processing framework. Stateful stream processing integrates the database and
the event-driven/reactive application or analytics logic into one tightly
integrated entity. With Flink, streaming data from many sources can be
ingested, processed, and distributed across various nodes. Flink can handle
graph processing, machine learning, and other complex event processing. Apache
Kafka888https://kafka.apache.org/documentation/streams/ is an open-source
publish and subscribe messaging solution. Services publishing (writing) events
to Kafka topics are asynchronously connected to other services consuming
(reading) events from Kafka - all in real-time. Kafka Streams lacks point-to-
point queues and falls short in terms of analytics. Spring Cloud Data
Flow999https://spring.io/projects/spring-cloud-dataflow is a microservice-
based streaming and batch processing platform. It provides tools to create
data pipelines for common use cases. Spring Cloud Data Flows has an intuitive
graphic editor that makes building data pipelines interactive for developers.
Amazon Kinesis Streams101010http://aws.amazon.com/kinesis/data-streams/ is a
service to collect, process, and analyse streaming data in real-time, designed
to get important information needed to make decisions on time. Cloud
Dataflow111111https://cloud.google.com/dataflow is a serverless processing
platform designed to execute data processing pipelines. It uses the Apache
Beam SDK for MapReduce operations and accuracy control for batch and streaming
data. Apache Pulsar is a cloud-native, distributed messaging and streaming
platform. Apache Pulsar121212https://pulsar.apache.org/ is a high-performance
cloud-native, distributed messaging and streaming platform that provides
server-to-server messaging and geo-replication of messages across clusters.
IBM Streams131313https://www.ibm.com/cloud/streaming-analytics proposes a
Streams Processing Language (SPL). It powers a Stream Analytics service that
allows to ingest and analyse millions of events per second. Queries can be
expressed to retrieve specific data and create filters to refine the data on
your dashboard to dive deeper.Source.

Discussion The real-time stream processing engines are developed for specific
use cases such as IoT, finance, advertisement, telecommunications, healthcare,
etc. Data is collected from some data-generating sources. Much data are of no
interest, and they can be filtered and compressed by orders of magnitude [8,
10, 11, 3]. They rely on distributed processing models, where unbounded data
streams are processed. Stream analytics is often performed after the complete
scanning of representative data sets. This strategy is inconvenient in real-
time stream processing to process the entire data stream at once. Windowing
mechanisms emerged for processing data stream in a predefined topology with a
fixed number of operations such as join, aggregate, filter, etc. The challenge
is to define these filters in such a way that they do not discard helpful
information.

Data streaming calls for research that can intelligently process raw data to a
size that its users can handle while not missing the needle in the haystack.
Furthermore, “online” analysis techniques must process streaming data on the
fly and combine them with historical data to provide past and current
analytics of observed environments.

##  3 H-STREAM for building querying pipelines for analysing streams

We propose H-STREAM, a Big stream operators platform implemented as micro-
services that can be composed for processing streams stemming from IoT
environments (see figure 1).

Figure 1: H-STREAM General Architecture.

H-STREAM is a middleware platform that provides Big Stream/Data operators
implementations like aggregation, descriptive statistics, filtering,
clustering, and visualisation as micro-services. Operators are composed to
define pipelines as queries that apply series of aggregation, filtering and
analytics operations to streams stemming from IoT farms and stored data. These
pipelines provide analytics visions of data to target applications. H-STREAM
relies on (i) message queues for collecting streams online from IoT farms; and
(ii) a backend execution environment [1] that provides a virtual data centre
with infrastructure resources necessary for executing costly processes.

Composing micro-services Micro-services can work alone or be composed to
implement simple or complex analytics pipelines (e.g., fetch, sliding window,
average, etc.). A series of data processing operations express queries in
H-STREAM applied to streams stemming from things, stores or micro-services. A
query is implemented by composing micro-services. The approach for composing
micro-services is based on a composition operation that connects them by
expressing a data flow (IN/OUT data). We currently compose aggregation
services (min, max, mean) with temporal windowing services (landmark, sliding)
that receive input data from storage support or a continuous data producer. We
propose connectors, namely Fetch and Sink micro-services that determine the
way micro-services exchange data from/to things, storage systems, or other
micro-services.

Figure 2: Taxonomy of queries that can be processed by composing micro-
services.

Big Stream Processing Pipelines Query Language We proposed a pseudo query
language with the syntax presented in figure 2, used to express:

  * •

The frequency in which data will be consumed (EVERY(number:Integer,
timeUnit:{minutes, seconds, hours})).

  * •

The aggregation function applies to an attribute of the input tuples (min,
max, mean).

  * •

The observation window on top of which aggregation functions will perform. The
window can involve only streams produced online (the last 5 seconds) or
include historical data (the last 120 days).

The observation can be done starting from a given instance in the past
(starting(number:Integer, timeUnit:{minutes, seconds, hours})) until something
happens. For example, from the instant in which the execution starts until the
consumer is disconnected. This corresponds to a landmark window [6]. It can
also be done continuously starting from a moving ”current instant” to several
{minutes, seconds, hours} before. This corresponds to a sliding window. The
expression includes the logic names of the data producers that can be a store
(Influx, Cassandra) and/or a stream queue provided by a message-oriented
middleware (e.g., RabbitMQ).

##  4 Stream processing micro-service

Figure 3 shows the general architecture of a stream micro-service. A micro-
service consists of three main components, Buffer Manager, Fetch and Sink, and
OperatorLogic. The micro-service logic is based on a scheduler that ensures
the recurrence rate in which the analytics operation implemented by the micro-
service is executed. Stream processing is based on “unlimited” consumption of
data ensured by the component Fetch that works if a producer notifies streams.
This specification is contained in the logic of the components OperatorLogic
and Fetch.

Figure 3: Architecture of a stream processing micro-service for processing
data streams.

As shown in the figure, a micro-service communicates asynchronously with other
micro-services using a message-oriented middleware. As data is produced, the
micro-service fetches and copies the data to an internal buffer. Then,
depending on its logic, it applies a processing algorithm and sends it to the
micro-services connected to it. The micro-services adopt the tuple oriented
data model as a stream exchange model among the IoT environment producing
streams and the micro-services. A stream is a series of attribute-value
couples where values are of atomic types (integer, string, char, float) from a
micro-service point of view. The general architecture of a micro-service is
specialised in concrete micro-services processing streams using well-known
window-based stream processing strategies: tumbling, sliding and landmark [7,
6]. Micro-services can also combine stream histories with continuous flows of
streams of the same type (the average number of connections to the Internet by
Bob of the last month until the next hour).

Since RAM assigned to a micro-service might be limited, and in consequence,
its buffer, every microservice implements a data management strategy by
collaborating with the communication middleware to exploit buffer space,
avoiding losing data and generating results on time. A micro-service
communicates asynchronously with other micro-services using a message-oriented
middleware. As data is produced, the micro-service fetches and copies the data
to an internal buffer. Then, depending on its logic, it applies a processing
algorithm and sends it to the micro-services connected to it. There are two
possibilities:

  * •

(i) on-line processing using tree window-based strategies [7, 6] (tumbling,
sliding and landmark) well known in the stream processing systems domain; (ii)

  * •

combine stream histories with continuous flows of streams of the same type
(the average number of connections to the Internet by Bob of the last month
until the next hour).

Figure 4 shows the general principle of the functional logic of a micro-
service dealing with windows that address the streams harvested before the
execution of the current query.

Figure 4: Synchronizing stream windows with historic data for computing
aggregations.

The following sections describe the strategies adopted for managing memory
coupled with cache and local disk to ensure the processing of every tuple
consumed by the micro-service; and how to compose micro-services for defining
data processing pipelines used by target applications for querying/analysing
data.

###  4.1 Interval oriented storage support for consuming streams

A micro-service that aggregates historical data and streams includes a
component named HistoricFetch. This component is responsible for performing a
one-shot query for retrieving stored data according to an input query (for
example, by a user or application). As described above, we have implemented a
general/abstract micro-service that contains a Fetch and Sink micro-services.
The historical fetch component has been specialized to interact with two
stores:

  * •

InfluxDB141414InfluxDB is a time series system accepting temporal queries,
useful for computing time tagged tuples (https://www.influxdata.comis) and

  * •

Cassandra151515Cassandra is a key-value store that provides non-temporal
read/write operations that might be interesting for storing huge quantities of
data (http://cassandra.apache.org).

The specification of the micro-service HistoricFetch and its specialisations
are briefly described in the following lines. The micro-service HistoricFetch
exports the following interface:
def queryToHistoric(function: String, value: String,
startTimeInMillis: Long, endTimeInMillis: Long,
groupByTimeNumber: Int, groupByTimeTimeUnit: String) :
List[List[Object]]

The method queryToHistoric(), shown in the code above, implements the
connection to a data store or DBMS, sends queries and retrieves data. It
returns the results packaged in the Scala structure List[ List[Row] ] objects,
where Row is a tuple of three elements:

  * •

Timestamp: Long, a timestamp in the format epoch;

  * •

Count: Double, the number of tuples (rows) that were grouped;

  * •

Result: Double, the result of the aggregation function.

As shown in the following code, a component HistoricFetch is created by the
micro-service specifying the name of the store (historicProvider), the name of
the database managing the stream history (dbName.series) and the execution
context. The current version of our micro-service runs on Spark, so the
execution context represents a Spark Context (sc).

      val hf : HistoricFetch = new
      HistoricFetch(historicProvider,dbName,series,sc)

The component HistoricFetch of the micro-service creates an object
HistoricProvider (step 1) that is used as proxy for interacting with a
specific store (i.e., InfluxDB or Cassandra). The store synchronously creates
an object Connection (step 2). The Connection object will remain open once the
query has been executed and results received by HistoricFetch. Then, the
component HistoricalFetch will use it for sending a temporal query using its
method queryToHistoric() (step 3). The result is then received in the variable
Result that is then processed (i.e., transformed to the internal structure of
the operator see below) and shared with the other components through the
Buffer of the micro-service (step 4).

Consider the observation of download and upload speed variations within users’
connections when they work on different networks. Assume that observations are
monitored online but that before the query is issued previous observations are
also stored. With a window and an average operators it is possible to answer a
query every two minutes give me the fastest download speed of the last 8
minutes. Figure 5 shows the corresponding expression in our language and it
shows how the logic of the operator deals with streams and data histories. The
query asks for a temporal interval filter that includes observations that
happened 8 minutes before the moment ”current” interpreted as the moment in
which the query starts its execution. It also asks to apply this temporal
filter every two seconds. This temporal filtering produces streams/data
batches on top of which the operator average is applied.

Figure 5: Big Streams Query Example.

###  4.2 Micro-services execution

Micro-services are executed on top of a Spark infrastructure deployed on a
virtual machine provided by the cloud provider Microsoft Azure (see Figure 6).
A micro-service running on top of the Spark platform processes streams
considering the following hypothesis:

  * •

There is a global time model for synchronizing different timelines (batch and
stream).

  * •

In the implementation, we use the Spark timeline as a global reference, and we
execute aggregations recurrently according to “time buckets”. The size of the
time bucket is determined by a query that defines an interval of observation
(e.g. the average number of connections of the last 5 hours) and a moving
temporal reference (e.g., every 5 minutes, until “now” if “now” is a moving
temporal reference).

Figure 6: Micro-service execution.

###  4.3 Experimental validation

We conducted experiments for validating the use of our micro-services. The
following lines describe the IoT environment that we prepared for our
experiments and some scale-up observations.

For deploying our demonstration experiment, we built an IoT farm using our
Azure Grant and implemented a distributed version of the IoT environment to
test a clustered version of Rabbit MQ and, therefore, address the scaling-up
problem in terms of the number of data producers (things) for our micro-
services.

Using Azure Virtual Machines (VM) we implemented a realistic scenario for
testing scalability in terms of:

  * •

(i) Initial MOM (RabbitMQ) installed in the VM2 in figure 7.

  * •

(ii) Producers (Things) installed in the VM1 in figure 7.

  * •

(iii) Micro-services installed in the VM3 in figure 7.

Figure 7: General experimental setting deployed on Windows Azure.

As shown in figure 7, in this experiment, micro-services and testbed were
running on separate VMs. This experiment leads to several cases scaling up to
several machines hosting until 800 things with a clustered version of Rabbit
using several nodes and queues that could consume millions of messages
produced at rates in the order of milliseconds see Figure 8).

Figure 8: Scale up scenarios

Observations in figure 9 showed the behaviour of the IoT environment regarding
the message-based communication middleware when the number of things
increased, when the production rate varies and when it uses one or several
queues for each consuming micro-service. We also observed the behaviour of the
IoT environment when several micro-services were consuming and processing the
data. Of course, the most agile behaviour is when the number of nodes and
virtual machine increases independently of the number of things. Indeed, note
that the performance of 800 things against 3 things does not change a lot by
increasing nodes, machines and queues. Note also that devoting one queue per
thing does not lead to essential changes in performance.

Figure 9: Scale up results

For our experimental scenarios, we varied the settings of the IoT environment
according to the properties characterising the scenarios. For the first one,
regarding the observation of neuroscientific experiments, we used fewer things
and queues and more nodes to achieve data processing in an agile way.
Regarding connectivity in cities with many people willing to connect devices
in different networks available in different urban spaces, we configured more
things and queues and nodes for the second one.

##  5 Use Cases

The objective of use cases was to have a proof of concept of the use of our
micro-services for implementing stream processing solutions to be weaved in
applications. They give insight into the way micro-services can be composed to
deal with streams production associated with volume and velocity. Then we
describe two use cases used for experimenting aspect (ii).

###  5.1 Experimental scenario I: Visual Analytics for Neuroscience

Neuroscience, like other experimental sciences, supports, refutes or validates
hypotheses by conducting experiments on living organisms. For instance, by
connecting electrical sensors to a cat’s spinal cord and monitoring its
neurons activities, neuroscientists can determine whether capsaicin (chilli
pepper active component) has the same effect as anaesthesia in the presence of
pain [9].

In a typical neuroscience experiment (see Figure 10), a neuroscientist is
responsible for: (i) preparing the subject; (ii) connecting and calibrating
sensors; (iii) collecting and storing the experiment data (e.g.file,
database); (iv) applying algorithms and statistics for discovering meaningful
patterns. Because of the volume of the resulting data and the complexity of
the algorithms used for finding patterns, the data analysis is usually done
post-mortem. Neuroscientists require novel tools for processing and exploring
data in real-time to better control the progress of an experiment. The reason
is that the experiment setting is complex (e.g., special and expensive
equipment, juridical protocols concerning experiments using animals, gathering
together field experts), long (e.g., 8 hours), and unique (e.g., every subject
has its proper characteristics).

Figure 10: Neuroscience scenario.

Therefore, we provided a solution for harvesting data produced during an
experiment and observing specific states of the vertebrae. Figure 11 gives an
overview of our approach. In the figure, data are collected during a
neuroscience experiment and continuously transmitted to our system. Data are
processed in real-time. Then, depending on the type of analysis that a
neuroscientist wants to conduct, she (i) defines queries using a set of
operators and (ii) chooses the kind of visual representation that she
requires. For instance, in our approach, a neuroscientist can group the data
into temporal windows of 1h. Then (for each window), she can choose different
visualizations (e.g., point chart, histogram, start plot) for analyzing the
correlation among the collected data.

Figure 11: Visual stream processing system for analyzing and visualizing data
streams.

We have conducted an experimental validation using data from studies regarding
pain performed in the neuroscience group at the Mexican research centre
CINVESTAV161616Special thanks to Diogenes Chavez from CINVESTAV Department of
Physiology, Biophysics and Neuroscience for providing the datasets used in
this work.. In this experiment, there is one micro-service for collecting the
data from an IoT environment (e.g. an a neuroscience experiment) and one for
plotting the data. The idea is that a consumer defines the sequence of micro-
services for processing and then plotting the data she desires to observe (as
shown in figure 12). For example, Give me the evolution of pain intensity in
L4ci or Give me the evolution of the minimum, average, maximum and intensity
of pain in L4ci every 3 seconds. Figure 12 shows how these queries are
implemented in terms of micro-services, including a window, aggregation, and
plotting ones. With this experiment, it was possible to observe online the
execution of the neuroscience experiment. Since streams were stored, it was
possible to observe the data from other experiments, compare them to make
decisions and adjust the phases of the experiment.

Figure 12: Observing continously the execution of the neurscience experiment.

###  5.2 Experimental scenario II: Analysing the behaviour of network services

The experimental scenario gives insight on the way micro-services can be
composed to deal with volume and velocity. For deploying our demonstration
experiment, we built an IoT farm using our Azure Grant and implemented a
distributed version of the IoT environment to test a clustered version of
RabbitMQ and, therefore, address the scaling-up problem in terms of several
data producers (things) for our micro-services.

The experimental scenario aims at analyzing the connectivity of the connected
society. The data set used for this demonstration has been produced in the
context of the Neubot project171717Neubot is a project devoted to measuring
Internet from the edges by the Nexa Center for Internet and Society at
Politecnico di Torino (https://www.neubot.org/).. It consists of network tests
(e.g., download/upload speed over HTTP) realized by different users in
different locations using an application that measures the network service
quality delivered by different Internet connection types 181818The HADAS group
previously used the Neubot data collection in the context of the FP7 project
S2EUNET in collaboration with Politecnico di Torino..

Figure 13: Analyzing the connectivity of connected society.

The idea is that people install the Neubot application on their computers and
devices. Every time they connect to the Internet using different networks (4G,
Ethernet, etc.), the application computes network quality metrics. as shown in
Figure 13.
The data is used then to answer queries such as:

  * •

The connection speed of one internet provider: Am I receiving the network
speed I am paying for all the time?

  * •

Internet availability with different Internet providers: Which are the periods
of the day in which I can upload/download files at the highest speed using
different network providers?

Figure 14: Micro-services architecture for analyzing society connectivity with
the neubot dataset.

The type of queries that we tested are the following. The result of an
H-STREAM query as a micro-services composition can be seen in Figure 15:

    EVERY 20 seconds compute the mean value of download_speed
          of the last 10 minutes
    FROM influxdb database neubot series speedtest and streaming
         RabbitMQ queue neubotspeed

Figure 15: H-STREAM query (micro-services composition) output.

    EVERY 60 seconds compute the max value of download_speed
          of the last 3 minutes
    FROM  cassandra database neubot series speedtests and streaming
          rabbitmq queue neubotspeed

    EVERY Ψ5 minutes compute the mean of the download_speed
            of the last 120 days
    FROM Ψcassandra database neubot series speedtests and streaming
            rabbitmq queue neubotspeed

    EVERY   30 seconds compute the mean value of upload_speed
            starting 10 days ago
    FROM Ψcassandra database neubot series speedtests and streaming
            rabbitmq queue neubotspeed

Depending on the observation window size, the micro-services access the
observations stored as post-mortem data sets and connect to online producers
currently observing their connections. For example, the third query wants to
observe a window of 10 days size. Our micro-services could deal with histories
produced in windows of size 10 days or even 120 days. Such massive histories
could be combined with recent streams and produce reasonable response times
(order of seconds).

##  6 Conclusions and Future Work

This paper proposes H-STREAM that provides micro-services for analysing and
exploring Big Streams in IoT environments. H-STREAM composes micro-services
deployed on high performance computing backends (e.g., cloud, HPC) to process
data produced by farms of things producing streams at different paces. Micro-
services compositions can correlate online produced streams with post-mortem
time series to discover and model phenomena or behaviour patterns of
intelligent environments. The advantage of this approach is that there is no
need for full-fledged data management services; micro-services compositions
can tailor data processing functions personalised to the requirements and the
characteristics of the applications and IoT environments.

Future work consists of developing a micro-services composition language that
can be used for expressing the data processing workflows that can be weaved
within target application logics [1]. We are particularly working on two urban
computing projects regarding the modelling and management of crowds and smart
energy management in urban clusters.

## References

  * [1] Akoglu, A., Vargas-Solar, G.: Putting data science pipelines on the edge. arXiv preprint arXiv:2103.07978 (2021)
  * [2] Alaasam, A.B., Radchenko, G., Tchernykh, A.: Stateful stream processing for digital twins: Microservice-based kafka stream dsl. In: 2019 International Multi-Conference on Engineering, Computer and Information Sciences (SIBIRCON). pp. 0804–0809. IEEE (2019)
  * [3] Apps, R., Scale, R.: Big Data Sourcebook (2011)
  * [4] Carbone, P., Ewen, S., Fóra, G., Haridi, S., Richter, S., Tzoumas, K.: State management in apache flink®: consistent stateful distributed stream processing. Proceedings of the VLDB Endowment 10(12), 1718–1729 (2017)
  * [5] Cardellini, V., Nardelli, M., Luzi, D.: Elastic stateful stream processing in storm. In: 2016 International Conference on High Performance Computing & Simulation (HPCS). pp. 583–590. IEEE (2016)
  * [6] Golab, L., Özsu, M.T.: Data stream management. Synthesis Lectures on Data Management 2(1), 1–73 (2010)
  * [7] Krämer, J., Seeger, B.: Semantics and implementation of continuous sliding window queries over data streams. ACM Transactions on Database Systems (TODS) 34(1), 4 (2009)
  * [8] Lyman, P., Varian, H.: How much information 2003? (2004)
  * [9] Martin, M., Contreras-Hernández, E., Béjar, J., Esposito, G., Chávez, D., Glusman, S., Cortés, U., Rudomin, P.: A machine learning methodology for the selection and classification of spontaneous spinal cord dorsum potentials allows disclosure of structured (non-random) changes in neuronal connectivity induced by nociceptive stimulation. Frontiers in neuroinformatics 9, 21 (2015)
  * [10] Woo, M.Y.: What’s the big deal about big data. Engineering and Science 76(3), 16–23 (2013)
  * [11] Zikopoulos, P., Eaton, C., et al.: Understanding big data: Analytics for enterprise class hadoop and streaming data. McGraw-Hill Osborne Media (2011)


# Protein folding is governed by memory-dependent friction

Benjamin A. Dalton  Fachbereich Physik, Freie Universität Berlin, 14195
Berlin, Germany  Cihan Ayaz  Fachbereich Physik, Freie Universität Berlin,
14195 Berlin, Germany  Lucas Tepper  Fachbereich Physik, Freie Universität
Berlin, 14195 Berlin, Germany  Roland R. Netz Corresponding author Fachbereich
Physik, Freie Universität Berlin, 14195 Berlin, Germany

###### Abstract

When described by a low-dimensional reaction coordinate, the rates of protein
folding are determined by free energy barriers, which separate folded and
unfolded states, and by friction, which results from interactions with the
solvent environment and internal, intra-molecular interactions. While it is
relatively straightforward to measure free energy barriers, a direct
evaluation of friction is far more elusive. Here, using memory-kernel
extraction methods, we directly calculate the friction acting on a protein,
which accounts for both external and internal contributions. By employing a
generalised Langevin equation (GLE) formalism to describe the dynamics of the
fraction of native contacts reaction coordinate $Q$, we evaluate friction
kernels for eight fast-folding proteins, taken from a published set of large-
scale, molecular dynamics protein simulations. Our results reveal that
friction is more important than the free energy barriers in determining
protein folding rates, particularly for larger proteins. We show that proteins
typically fold in an accelerated kinetics regime, which is a purely non-
Markovian effect.

## Introduction

In many cases, protein functionality depends on proteins folding successfully
into a specific three-dimensional conformational state. Protein folding
requires that a linear polypeptide chain is driven to explore conformation
space by interactions with a solvating environment, and is shaped by both
solvent interactions and internal interactions between amino acids. Folding
and unfolding typically occur over a free energy landscape [1, 2, 3, 4, 5, 6],
with distinct states separated by free energy barriers, often of the order of
a few $k_{\rm{B}}T$’s in height. Theories for describing protein folding are
often phrased as reaction rate theories with an explicit dependence on the
free energy barrier of separation, and friction. An accurate understanding of
the kinetics of barrier-crossing reactions is not only important for
describing protein folding but is essential in many other fields, such as
nucleation theory and chemical kinetics. Reaction rate theory dates back to
Arrhenius, who, in 1889, [7] showed that the transition times between reactant
and product scale according to $\kappa\text{e}^{U_{0}/k_{\rm{B}}T}$, where
$U_{0}$ is the height of the energy barrier separating the two states. The
pre-factor $\kappa$, however, depends on both the free energy profile and
friction. Building on the transition state theory of Eyring, [8], Kramers was
the first to include an explicit dependence on solvent friction [9]. The
friction in Kramers theory is frequency-independent, suggesting that the
environment of the reacting system relaxes infinitely fast compared to the
relevant time scales of the reacting system. Regardless of its simplicity,
Kramers theory predicts reaction rates in both the over-damped and inertia-
dominated regimes well and is sufficient for describing many systems. Various
advancements have been made that bridge the over-damped and inertia-dominated
regimes, with many accommodating the finite environmental relaxation times
[10, 11, 12, 13, 14]. In the case of proteins, it is known that friction is
determined by a combination of interactions with the solvent environment and
internal intramolecular interactions [15, 16, 17], and that finite relaxation
is occurs [18]. A direct evaluation of the friction acting on a protein has
not been possible and approaches have rather relied on determining friction
from memoryless reaction rate theory [19, 20, 21]. In this paper, we directly
extract the friction acting on the conformational dynamics of proteins from
molecular dynamics simulation trajectories, which accounts for both internal
and external friction contributions, and finite relaxation times.

We evaluate friction for eight fast-folding proteins, obtained from previously
published large-scale molecular dynamics simulation trajectories [22]. These
cutting-edge simulations, performed by the Shaw group using the purpose-built
Anton super-computer [23, 24], represent a breakthrough in the simulation-of-
scale, enabling all-atom simulation of proteins that would not otherwise be
practically possible. The lengths of the proteins range between 10 and 80
amino-acid residues, and the simulation times between approximately 100 $\mu$s
and 3 ms. The result is a diverse protein data set, comprised of mixtures of
$\alpha$-helix and $\beta$-hairpin secondary structures, and an assortment of
tertiary structures, where all proteins execute a multitude of folding and
unfolding events.

We analyse all protein trajectories in the framework of a non-Markovian,
generalised Langevin equation (GLE). The GLE is a low-dimensional
representation of some higher dimension system. In the present case, we
collapse the all-atom dynamics of the composite water-protein systems onto a
one-dimensional reaction coordinate. This process of dimension reduction is
known as projection [25, 26]. We project the all-atom trajectories, as
provided by the Shaw group, onto a well-known fraction of native contacts
reaction coordinate [27] with soft cut-off $Q$ [28], and hence extract
friction memory kernels from $Q(t)$. Free-energy profiles and friction are a
property of dimension reduction, and are therefore unique to a given reaction
coordinate. The extracted friction kernels encode dissipation across a spread
of time scales, representing the finite relaxation rates of solvent processes
and internal reconfigurations, as far as they are represented in $Q(t)$.
Consequently, the friction kernel contains the information for the full
friction acting on the reaction coordinate. From the free energy landscape and
friction kernels, we derive predictions for the folding times on different
levels of reaction rate theory and compare to the folding times measured in
the simulation. In doing so, we show that the best prediction for the measured
folding times is given by a multi-time-scale, non-Markovian theory, which has
hitherto only been applied to the characterisation of model systems [13, 14].
This validation is only possible owing to the long simulation times and
diverse range of proteins that comprise the extensive protein data set. We
show that the memory decay-times, i.e., the duration of memory effects, are
significantly long for all proteins in this set, in some instances as long as
the folding times. This indicates that even for $Q$, which is typically
considered by other measures to be a good reaction coordinate [28], must be
considered as a poor reaction coordinate, when judged according to its non-
Markovianity. Finally, we also show that, for this particular set of proteins,
friction is more important than free-energy barriers in determining folding
times, and that this dominance of friction increases for larger proteins.
Taken together, our findings suggest that, when represented by a low-
dimensional reaction coordinate, the conformational dynamics of proteins are
dominated by friction and should be treated as non-Markovian in nature.

## Results

Figure 1: The folding and unfolding of eight fast-folding proteins. A) PDB
entries for eight proteins, with the number of amino acids in each protein. B)
250 $\mu$s trajectory segment for the $Q(t)$ reaction coordinate
($\alpha_{3}$D protein). Left magnification: a sequence of folding first-
passage times, from the unfolded state $Q_{\rm{u}}$ to the barrier top
$Q_{\rm{b}}$. Right magnification: an example folding transition path and
corresponding transition path time. C) Free energy profile for the
$\alpha_{3}$D protein. Configuration snapshots show example folded and
unfolded states, with an example barrier state. $Q_{\rm{u}}$ and $Q_{\rm{b}}$
are the reaction coordinate values in the folded and unfolded states,
respectively. $Q_{\rm{b}}$ is the transition barrier top. Asymmetric barrier
heights are also indicated. For $\alpha_{3}$D, the barrier faced by the
unfolded protein,
$U^{\rm{u}}_{0}=U(Q_{\rm{b}})-U(Q_{\rm{u}})=3.2\;k_{\text{B}}T$, is less than
the barrier faced by the folded protein
$U^{\rm{f}}_{0}=U(Q_{\rm{b}})-U(Q_{\rm{f}})=1.8\;k_{\text{B}}T$. The distance
from the unfolded state to the barrier top $L_{\rm{u}}=Q_{\rm{b}}-Q_{\rm{u}}$
is greater than the distance from the folded state to the barrier top
$L_{\rm{f}}=Q_{\rm{f}}-Q_{\rm{b}}$. D) Normalised GLE memory kernel, extracted
from $Q(t)$ for the $\alpha_{3}$D protein. Running integral $G(t)$ and the
limiting total friction $\gamma$ (inset - dashed line). Running integral of an
exponential series fit to $\Gamma(t)$ (red curve). E) Total friction $\gamma$
for each protein, plotted as a function of the number of residues $N$, divided
by $k_{\text{B}}T$ to account for unique temperature $T$ (power law with
exponent $\sim 2.8$ (red line)). F) Effective mass $m$, plotted as a function
of $N$ (linear fit with a gradient of 41 u per residue (red line)). G) and H)
show the Arrhenius factor $\text{e}^{\text{U}_{0}/k_{\text{B}}T}$ and barrier
crossing times $\tau_{\rm{MFP}}^{\text{MD}}$ as a function of protein lengths
folding and unfolding transitions individually.  Figure 2: Comparison of
relevant time scales for protein folding and unfolding dynamics and kinetics.
A) The inertia time scales $\tau_{\rm{m}}$ compared to the diffusion times
$\tau_{\rm{D}}$, showing that all systems are in the over-damped regime. B)
Memory decay time-scale $\tau_{\rm{mem}}$, calculated from the first moment of
the memory kernel $\Gamma(t)$, compared to the diffusion times $\tau_{\rm{D}}$
. The light grey lines indicate the bounding regime for $\tau_{\rm{mem}}$
between $\tau_{\rm{D}}{\times}10^{-2}$ and $\tau_{\rm{D}}{\times}10^{1}$,
which is the domain in which memory-induced kinetic speed-up is expected. C)
Memory times $\tau_{\rm{mem}}$ compared to the folding and unfolding times,
expressed as the mean first passage times $\tau_{\rm{MFP}}^{\text{MD}}$. D)
Memory times compared to the transition path times $\tau_{\rm{MTP}}$ leading
from the folded and unfolded state minima to the barrier tops. The broken
lines in each plot indicate exact equivalence, or correlation, between the
respective quantifies.

Extracting friction from protein simulations: The PDB structures for the eight
proteins, along with the number of residues that make up each chain, are shown
in Fig. 1A. The fraction of native contacts reaction coordinate $Q$ is
described in detail in the Materials and Methods section. In brief, $Q$
provides a measure for the deviation of a given state away from some native
reference state, which typically represents the folded state. $Q$, which is
evaluated in connectivity space and is therefore unitless, is closer to $Q=1$
in the folded state and closer to $Q=0$ when unfolded. For the set of eight
proteins, we identify a native state and hence project each corresponding all-
atom trajectory onto a unique $Q(t)$. In Fig. 1B, we show a 250 $\mu$s
trajectory segment for the $\alpha_{3}$D protein. The distinct folded and
unfolded states are discernible, located at $Q_{\rm{f}}=0.75$ and
$Q_{\rm{u}}=0.51$ respectively, separated by a barrier at $Q_{\rm{b}}=0.66$.
In the magnification windows set above the trajectory, we indicate schematic
examples for the calculation of first-passage times $\tau_{\rm{FP}}$ and
transition path times $\tau_{\rm{TP}}$. Throughout this paper, we consider
folding transitions as leading from the unfolded state minimum to the barrier
top ($Q_{\rm{u}}\rightarrow Q_{\rm{b}}$), and unfolding transitions as leading
from the folded state minimum to the barrier top ($Q_{\rm{f}}\rightarrow
Q_{\rm{b}}$). This accounts for the strong asymmetries observed in the free
energy profiles for the eight proteins. In the left magnification, we identify
a sequence of folding first-passage events. The mean of all such first-passage
events, taken over the full trajectory, are accumulated to evaluate the mean
first-passage time $\tau_{\rm{MFP}}^{\text{MD}}$ for the folding reaction.
Likewise, we evaluate the unfolding times as the mean of all first-passage
events from $Q_{\rm{f}}$ to $Q_{\rm{b}}$. In the right magnification, we show
an example folding transition path, connecting the unfolded state minimum to
the barrier top. A mean transition path time $\tau_{\rm{MTP}}$ is simply the
mean for all such transition path durations for either folding or unfolding
transition. From the trajectories of $Q(t)$ for each protein, we calculate the
free energy profile $U(Q)$ acting on the reaction coordinate. In Fig. 1C, we
show $U(Q)=-k_{\rm{B}}T\text{log}[\rho(Q)]$ for the $\alpha_{3}$D protein,
where $\rho(Q)$ is the probability density for $Q(t)$, $k_{\rm{B}}$ is
Boltzmann’s constant, and $T$ is the system temperature, which has a unique
value for each protein. We show the free energy profiles for all proteins in
the SI (Section 6). All free energy profiles are asymmetric, such that the
barrier heights faced by the folded and unfolded states ($U_{0,\rm{u}}$ and
$U_{0,\rm{f}}$) are not equal. Likewise, the distances in reaction-coordinate
space from the minima to the barrier top ($L_{\rm{f}}$ and $L_{\rm{u}}$) are
also different.

Central to this paper is the extraction of time-dependent friction kernels
$\Gamma(t)$ from the $Q(t)$ trajectories for each of the eight proteins. We
describe the stochastic dynamics of $Q(t)$ by the one-dimensional generalised
Langevin equation (GLE):

$\begin{split}m\ddot{Q}(t)&=-\int\limits_{0}^{t}\Gamma(t-t^{\prime})\dot{Q}(t^{\prime})dt^{\prime}\\\
&\quad\quad\quad\qquad-\nabla U\big{[}Q(t)\big{]}+F_{R}(t),\end{split}$ (1)

where $\Gamma(t)$ is the friction memory kernel, $F_{R}(t)$ is the random
force term satisfying the fluctuation-dissipation theorem $\langle
F_{R}(t)F_{R}(t^{\prime})\rangle=k_{\rm{B}}T\Gamma(t-t^{\prime})$, $m$ is the
effective mass of the reaction coordinate, and $U(Q)$ is the free energy
profile. Eq. 1 neglects non-linear friction and is therefore approximate [29].
There are various methods for extracting dynamic friction from discrete time-
series trajectories [30, 31, 32]. We describe the method used here for
extracting $\Gamma(t)$ from $Q(t)$ in detail in the SI (Section 2). In short,
we directly extract the running integral of the memory kernel
$G(t)=\int_{0}^{t}\Gamma(t^{\prime})dt^{\prime}$ using a Volterra extraction
scheme, which is suitable for arbitrary, non-linear, free-energy profiles
$U(Q)$ [33, 34]. The memory kernel is then given by $\Gamma(t)=dG(t)/dt$,
which is evaluated numerically. The memory kernel for the $\alpha_{3}$D
protein, normalised by $\Gamma(0)$, is shown in Fig. 1D. The dynamics of the
$\alpha_{3}$D protein exhibit significant memory effects, which is clear from
the inset, which shows that the running integral $G(t)$ plateaus after about
10 $\mu$s. We show the memory kernels for all proteins in the SI (Section 7).

Scaling of GLE parameters with system size: Having extracted the full time-
dependent friction kernel via $G(t)$, we also have access to the zero-
frequency friction $\gamma=G(t\rightarrow\infty)$, which we refer to
throughout as the total friction. In Fig. 1E, we show $\gamma$ for each
protein, plotted as a function of the number of residues in each chain $N$. To
account for each system having a unique temperature, we divide by
$k_{\text{B}}T$. We fit a power law and find that
$\gamma/k_{\text{B}}T=3.2{\times}10^{-3}N^{2.8}\;\mu\text{s}$ (red line) [35].
If the friction on $Q$ is determined solely by solvent interactions, we would
expect approximately linear scaling. This is the case for the effective mass,
which we calculate using the equipartition theorem
$m=2k_{\text{B}}T/\langle\dot{Q}\rangle^{2}$ (Fig. 1F). The linear scaling of
$m$ (red line) indicates that the effective mass of $Q$ is largely determined
by the molecular mass of the protein chain [29]. The power-law scaling of
$\gamma$ suggests that intramolecular interactions are more influential for
longer chains, owing to the reptation-like reconfiguration dynamics [36].
Interestingly, increasing chain length does not appear to directly couple to
free energy barrier heights. The Arrhenius barrier terms
$\text{e}^{\text{U}_{0}/k_{\text{B}}T}$, plotted in Fig. 1G, show that the
hindrance to protein folding imposed by free energy barriers is remarkably
uncorrelated with the size of a protein. The folding and folding kinetics,
however, which we quantify as the mean first-passage times
$\tau_{\rm{MFP}}^{\rm{MD}}$ (Fig. 1H), do increase with $N$. Therefore,
whereas the separation barrier heights do not appear to increase with $N$, the
barrier crossing times do. Overall, Fig. 1E-H tells us that, when projected
onto $Q$, friction plays a dominant role in governing the dynamics and
kinetics of protein folding and that the free energy profiles, while certainly
also essential, appear to be less influential than friction in determining
protein folding reaction rates, as will be correlated more quantitatively
further below.

Figure 3: Predicting folding and unfolding times. A) Diffusion times
$\tau_{\rm{D}}$, according to Eq. 2, compared to the measured barrier crossing
times $\tau_{\rm{MFP}}^{\rm{MD}}$. B) Markovian predictions
$\tau_{\rm{MFP}}^{\rm{Mar}}$, according to Eq. 3, for barrier crossing times
on exact free-energy profiles. C) Free energy factor $\xi_{\text{U}}$,
accounting for the effects of the free-energy landscape on the reaction time,
according to Eq. 3. D) Non-Markovian predictions
$\tau_{\rm{MFP}}^{\rm{noMar}}$, according to Eq. 5, for the barrier crossing
times, evaluated by correcting $\tau_{\rm{MFP}}^{\rm{Mar}}$ to account for
non-Markovian effects, using the correction factor $\xi_{\text{noMar}}$, as
given by Eq. 4. Root-mean-squared logarithmic deviations score the population
deviation from exact agreement with the measured $\tau_{\rm{MFP}}^{\rm{MD}}$
for A), B), and D) (dashed lines), and the optimal linear regression in C)
(red line). See Fig. 2 for symbol legend.

Memory duration is significant in protein folding: As well as proving the
total friction for $Q$, the extraction of $\Gamma(t)$ provides a time scale
for the sustain of memory effects. Furthermore, $\gamma$ can be used
predictively, to evaluate other key dynamic time scales. How do these time
scales compare? The inertia time $\tau_{\text{m}}=m/\gamma$ is the time scale
beyond which the system becomes diffusive. The diffusion time:

$\tau_{\rm{D}}=\frac{\gamma L^{2}}{k_{B}T},$ (2)

is the time taken for a Brownian particle to diffuse over a characteristic
distance $L$ in the absence of free energy gradients. Both $\tau_{\text{m}}$
and $\tau_{\rm{D}}$ depend on $\gamma$. A system is in an over-damped regime
when $\tau_{\text{m}}\ll\tau_{\rm{D}}$. This condition is met for all eight
proteins, as can be seen in Fig. 2A, for both $L_{\rm{u}}$ and $L_{\rm{f}}$.
One ambiguous case is Chignolin unfolding, for which
$\tau_{\rm{m}}/\tau_{\rm{D}}=0.027$. For Chignolin, both the characteristic
length in the folded state $L_{\rm{f}}$, and $\gamma$, are considerably small.

By comparing $\tau_{\rm{D}}$ to the decay times of the extracted memory
kernels, we can intuit whether we expect memory effects to influence barrier
crossing processes. Memory-dependent friction typically exhibits cascading
time scales, meaning that $\Gamma(t)$ is multi-modal, with time scales
spanning many orders of magnitude [33]. To assign a single effective time-
scale to the memory decay $\tau_{\rm{mem}}$, we evaluate the first-moment for
the set of extracted memory kernels:
$\tau_{\rm{mem}}=\int_{0}^{\infty}t\Gamma(t)dt/\int_{0}^{\infty}\Gamma(t)dt$.
In the SI (Section 3), we address the issue of discretization due to the
intrinsically low time resolution of the MD data with a highly resolved
trajectory for an $\alpha$-helix forming Ala9 homo-peptide chain. Fig. 2B
shows that overall, memory times are non-negligible when compared to diffusion
times. Memory times accelerate barrier crossing kinetics in the range
$1{\times}10^{-2}<\tau_{\rm{mem}}/\tau_{\rm{D}}<1{\times}10^{1}$ [12, 13, 14],
which is indicated on Fig. 2B, where the Grote-Hynes theory is valid [11]. In
Fig. 2C, we see that memory times are at most comparable to the reaction
times, but are typically shorter, indicating that $Q$ is a poor reaction
coordinate. The overall concurrence between $\tau_{\rm{mem}}$,
$\tau_{\rm{D}}$, and $\tau_{\rm{MFP}}^{\text{MD}}$ reveals a coupling between
memory effects and folding kinetics.

The comparison between memory times and transition path times is particularly
revealing. Transition paths are the segments of a trajectory where the protein
actually executes the reconfiguration from one specific state to another
target state. Transition paths are of much interest since they provide
information about the actual folding mechanisms of a protein, and have been
studied extensively in experiments [37, 38, 39, 40], and in the context of
simulation and theory [22, 28, 41, 42, 43, 44, 45]. In Fig. 2D, we show that
for all proteins, memory times are long compared to transition path times.
This is significant for two reasons. Firstly, it tells us that entire
transitions either fold or unfold under the influence of memory effects that
linger from the preceding state. Secondly, it explains why the approximate,
position-independent form of Eq. 1 works. Even describing the dynamics of $Q$
as having distinct memory dependence on separate sides of the barrier would be
neutralised since memory effects can last for as much as two orders of
magnitude longer than the time required to move between states. Therefore, we
must keep in mind the overwhelming dominance of the duration of memory effects
compared to the time required to traverse a free energy barrier.

Figure 4: Barrier crossing times indicate memory-induced speed-up. A)
Deviations from Markovian barrier crossing kinetics, plotted as a function of
scaled memory times $\tau_{\rm{mem}}/\tau_{\rm{D}}$, for the folding and
unfolding of all eight proteins. Rescaled MD simulation values for each
protein ($\tau_{\rm{MFP}}^{\text{MD}}/\tau_{\rm{MFP}}^{\rm{Mar}}$, red and
black symbols) are compared to multi-modal, non-Markovian prediction
($\xi_{\text{noMar}}$, blue and yellow symbols). See Fig. 2 for symbol legend.
B) Multi-modal, non-Markovian predictions $\xi_{\text{noMar}}$ plotted as a
function of rescaled first-moment memory time $\alpha\tau_{\rm{mem}}$ for two
examples: Trp-Cage and $\alpha_{3}$D. The symbols are for $\alpha=1$, which
coincide with the appropriate $\xi_{\text{noMar}}$ values in A). Dashed lines
indicate quadratic scaling
$\sim\big{(}\tau_{\rm{mem}}/\tau_{\rm{D}}\big{)}^{2}$ for the long memory-time
slow-down regime.

Non-Markovian models predict protein folding times: From the simulation
trajectories, we measure the folding and unfolding times, evaluated via the
mean first-passage times $\tau_{\rm{MFP}}^{\text{MD}}$. Since we can also
extract the total friction directly from the simulation trajectories, we can
parametrise various theoretical estimates for folding and unfolding times, and
hence compare to the measured simulation values. To quantify the population-
wide deviation away from the measured $\tau_{\rm{MFP}}^{\text{MD}}$, we use
the root-mean-square logarithmic deviation (RMSLD). We collect the folding and
unfolding times into a single set and we logarithmically weight each value to
balance contributions from the spread of time scales (see SI - Section 10).
The population-wide deviations are given by
$\text{RMSLD}=\sqrt{\sum_{i=1}^{N}(\text{log}(\tau^{\text{MD}}_{\text{MFP},i})-\text{log}(\tau^{\text{theo}}_{\text{MFP},i}))^{2}/N}$,
where $\tau^{\text{MD}}_{\text{MFP},i}$ is the measured value from MD,
$\tau^{\text{theo}}_{\text{MFP},i}$ is the theoretical predicted value for a
given method, and $N=16$ is the total number of reactions. For the diffusion
times $\tau_{\rm{D}}$, we determine a total logarithmic deviation of 1.51.
This result is remarkable since it completely neglects the free energy
barriers, depending only on the friction and characteristic lengths. This
result implies that friction alone provides a good approximation for protein
folding times. A Markovian prediction, which explicitly accounts for the
extracted free energy profiles, is given by:

$\begin{split}&\tau_{\rm{MFP}}^{\rm{Mar}}\big{(}Q_{\text{a}}|Q_{\text{b}}\big{)}\\\
&\qquad=\tau_{\rm{D}}\int_{Q_{\text{a}}}^{Q_{\text{b}}}\text{e}^{\beta
U(x)}\bigg{[}\int_{-\infty}^{x}\text{e}^{-\beta
U(y)}\frac{dy}{L}\bigg{]}\frac{dx}{L}\\\
&\qquad=\tau_{\text{D}}\xi_{\text{U}},\end{split}$ (3)

where $\beta=1/k_{\rm{B}}T$. Eq. 3 assumes constant friction across $Q$ and
evaluates the mean first-passage times between any two points on the free
energy landscape, $Q_{\text{a}}$ and $Q_{\text{b}}$. In the case of folding
reactions, we calculate
$\tau_{\rm{MFP}}^{\rm{Mar}}(Q_{\text{u}}|Q_{\text{b}})$. For unfolding, we
reflect the free energy profile $U(Q)\rightarrow U(Q_{\rm{max}}-Q)$ and hence
calculate $\tau_{\rm{MFP}}^{\rm{Mar}}(Q_{\text{f}}|Q_{\text{b}})$. Again, we
see that the population-wide prediction is good (Fig 3B), but with a general
trend that the Markovian predictions are slow compared to the MD data. From
Eq. 3, $\xi_{\text{U}}$ is a factor that depends on the exact free energy
profile. The RMSLD value of 1.41 for the Markovian predictions shows that the
contribution due to $\xi_{\text{U}}$ is small. We emphasised this in Fig 3C,
where we show the population of $\xi_{\text{U}}$ values corresponding to the
set of barrier crossing times. Here, we calculate an RMSLD value of 1.98 using
linear regression (red line). This large spread of $\xi_{\text{U}}$ values
confirms that by explicitly including free energy profiles, we can not improve
upon a prediction for the barrier crossing times that is given solely by
friction.

Can we improve these predictions by taking the proposed non-Markovianity into
account? Using a simplified non-Markovian model with multi-modal memory in
symmetric free energy profiles with parabolic barriers [13, 14], we determine
a non-Markovian correction factor $\xi_{\text{noMar}}$, which is given by:

$\xi_{\text{noMar}}=\frac{\tau^{\text{H}}_{\rm{MFP}}(\tau_{\rm{D}},U_{0},\\{\gamma_{i}\\},\\{\tau_{i}\\},\tau_{\rm{m}})}{\tau_{\rm{MFP}}^{\text{H}}(\tau_{\rm{D}},U_{0},\\{\gamma_{i}\\},0,0)}.$
(4)

$\tau^{\text{H}}_{\rm{MFP}}$ is the barrier-crossing time predicted by a
heuristic cross-over formula. This term contains a series of $M$ over-damped
contributions, which dominate when memory times are short compared to the
diffusion times, and $M$ energy-diffusion contributions, which dominate when
memory times are large. The exact nature of these contributions is described
in detail in the materials and methods section (see Eqs. M2 \- M4). Briefly,
we fit the extracted memory kernels with an exponential series
$\Gamma(t)\approx\sum_{i=1}^{M}(\gamma_{i}/\tau_{i})\exp{(t/\tau_{i})}$ (see
SI - Section 7), resulting in a set of $M$ amplitudes $\\{\gamma_{i}\\}$ and
time-scales $\\{\tau_{i}\\}$. For Chignolin, $M$=2. For all other proteins,
$M$=3. The denominator in Eq. 4 represents the memoryless, Markovian limit,
where the mass is assumed to be zero and the memory time scales are infinitely
short. The results for the corrected Markovian predictions:

$\tau_{\rm{MFP}}^{\text{noMar}}=\xi_{\text{noMar}}\tau_{\rm{MFP}}^{\rm{Mar}},$
(5)

are given in Fig 3D. We see an overall improvement, except for the case of
Chignolin unfolding, Villin folding, and Protein G unfolding. Population-wide,
we see a significant improvement in the RMSLD, yielding a value of 0.86. While
the errors are still quite large, we recall that aside from the sets
$\\{\gamma_{i}\\}$ and $\\{\tau_{i}\\}$, which are taken from the extracted
memory kernels, $\xi_{\text{noMar}}$ is a fit-free quantity, with all
parameters extracted from the MD trajectories. Non-Markovian reaction rate
theory represents a significant improvement when predicting protein folding
times.

Proteins fold in the memory-induced speed-up regime: In reactive systems,
memory effects can result in either the speed-up or slow-down of reaction
times, when compared to a memoryless analog, depending on the memory time
scales [12, 13, 14]. For our analysis, $\tau_{\rm{MFP}}^{\rm{Mar}}$ represents
a memoryless analogue to the actual, measured system
$\tau_{\rm{MFP}}^{\text{MD}}$. The ratio of
$\tau_{\rm{MFP}}^{\text{MD}}/\tau_{\rm{MFP}}^{\rm{Mar}}$ reveals deviations
from pure Markovian behaviour. If we rescale memory times by the corresponding
diffusion times, we have a relative measure for the duration of memory decay
compared to the underlying transport time over the reactive distance. By
plotting $\tau_{\rm{MFP}}^{\text{MD}}/\tau_{\rm{MFP}}^{\rm{Mar}}$ against the
rescaled memory times $\tau_{\rm{mem}}/\tau_{\rm{D}}$ for all proteins, we
observe the deviation from pure Markovianity over the range of extracted
memory times. The red (folding) and black (unfolding) symbols in Fig. 4A
reveal a speed-up regime across the population of proteins. In the short and
long memory-time limits, reaction times trend to Markovian behaviour, which is
expected from previous works [13, 14]. Memory-induced speed-up has already
been observed in short homo-peptide chains [33]. Our results indicate that
memory effects are present in larger proteins and that these effects likely
explain the observed acceleration of the folding and unfolding reactions.

Using $\xi_{\text{noMar}}$ predictively, we can demonstrate that the reaction-
acceleration effects are induced by memory effects. $\xi_{\text{noMar}}$, as
given by Eq. 4, is an estimate for a reaction time, predicted by a simple non-
Markovian model, normalised by the over-damped, Markovian limit. We plot the
parametrised predictions $\xi_{\text{noMar}}$ for the set of proteins over the
rescaled memory times (blue (folding) and yellow (unfolding) symbols in Fig.
4A). Except for the case of Chignolin unfolding, which is predicted to be much
slower than the measured value, the agreement between the predicted folding
times and those measured in the simulations is good, indicating that a simple,
multi-modal, non-Markovian model describes well the accelerated barrier-
crossing observed in protein folding simulations.

Reaction speed-up effects are predicted by the Grote-Hynes theory, which does
not account for slow-down effects as $\tau_{\rm{mem}}\rightarrow\infty$, known
to occur in model systems. We can contextualise the data for $\xi$, as given
in Fig. 4A, by uniformly rescaling all time scales that enter into the
heuristic formula Eq. 4 by a factor $\alpha$, while keeping all other
parameters fixed. Such a rescaling might reflect a variation in solvent
viscosity, for example. In doing so, we effectively generate an
$\alpha$-dependent memory kernel
$\Gamma^{\alpha}(t)=\sum_{i=1}^{M}(\gamma_{i}/\alpha\tau_{i})\exp{(t/\alpha\tau_{i})}$,
and a corresponding effective single-valued time-scale
$\alpha\tau_{\rm{mem}}=\int_{0}^{\infty}t\Gamma^{\alpha}(t)dt/\int_{0}^{\infty}\Gamma^{\alpha}(t)dt$.
In Fig. 4B, we show two examples of $\xi(\alpha)$, plotted as a function of
$\alpha\tau_{\rm{mem}}/\tau_{\rm{D}}$. The example of Trp-Cage is particularly
interesting as it shows that, in the context of this rescaling curve, the
prediction for unfolding, as given in Fig. 4A and replotted as the black
circle in Fig. 4B, actually coincides with the inflection from the speed-up
regime, i.e. $\xi_{\text{noMar}}<1$, to the kinetic slow-down regime,
$\xi_{\text{noMar}}>1$. This is also true of Villin and the WW-Domain (see SI
- Section 8). The kinetic slow-down regime is characterised by quadratic
scaling in the long memory-time limit, which is also displayed for the
$\alpha$ rescaling. Overall, the rescaling curves show that all proteins
folded and unfold in the kinetic speed-up regime, which is a purely non-
Markovian effect.

## Discussion

We have extracted memory kernels from large-scale MD simulation trajectories
of eight fast-folding proteins and found that, when measured in the $Q$
reaction coordinate, memory times are comparable to the protein folding and
unfolding reaction times, and typically vastly exceed the transition path
times. It is unclear what makes a good reaction coordinate [46, 47], and such
a quality is often context-dependent. For example, a reaction coordinate might
be optimized such that barrier-crossing states are most likely to be
transition states [48]. $Q$ has already been shown to exhibit this property
[28]. But such definitions do not say anything about dynamics. In the SI
(Section 9), we show examples of $\Gamma(t)$ for other standard reaction
coordinates, such as the end-to-end distance, radius of gyration, and the RMSD
from the native state. The resultant memory times span orders of magnitudes
but are typical of the order of the diffusion time scale $\tau_{\text{D}}$.
Therefore, none of these reaction coordinates can be considered a good
reaction coordinate, insofar as optimizing for Markovianity is a viable
measure of goodness. However, we have shown that it does not matter whether we
have a good or bad reaction coordinate. So long as we have the appropriate
non-Markovian framework for analysing protein folding trajectories, taken from
either simulation or experiments [37, 38, 40, 49, 50], we can accurately
predict the folding kinetics of a protein. A simple multi-component, non-
Markovian model is the most accurate for predicting the folding time’s measure
in the MD simulations.

Our results indicate that, for the set of proteins considered in this paper,
friction is more important than free energy barriers in determining folding
times. The disclaimer is important here since the protein trajectories are
drawn from a biased set, such that the proteins are selected for their short
folding times, and simulation temperatures are chosen to maximize
folding/unfolding transitions, i.e., the melting temperatures. Our results
show that the total friction acting on $Q$ increases more than quadratically
with chain length. Barrier heights do not increase with chain length, but the
barrier crossing time do (Fig. 1E-H). We cannot suggest that such dependencies
hold for proteins with much larger barrier heights, or for general
temperatures. Regardless, it appears that, for the present data, contributions
from the friction-dependent pre-factor in the Arrhenius-type reaction-rate
theories for protein folding dominate the exponential term. This is made clear
when we compare Figs. 3A - C. A simple Markovian prediction for the barrier
crossing times, one that only depends on friction, i.e. $\tau_{\text{D}}$,
represents well the measure barrier crossing times. Explicitly accounting the
exact free energy profiles makes very little improvement on this prediction.
The friction is dominant. A far greater improvement is made when we account
for the non-Markovian effects, as is seen in Fig. 3D.

In Eq. 3, we have only assumed constant friction. Alternative models include
position-dependent friction. In the SI (Section 4), we show that for a
Markovian model with position-dependent friction, there does not exist a
unique friction profile $\gamma(Q)$ for describing both folding and unfolding
kinetics. This was already shown to be true for a simple $\alpha$-helix
forming homo-alanine chain. From this we conclude that we can not predict both
folding and unfolding times consistently for a given $\gamma(Q)$.

Whether in simulation or experiment, protein conformation dynamics must be, in
general, described by some reduced-dimension collective reaction coordinate.
Overall, our results suggest that, irrespective of the choice of the reaction
coordinate, non-Markovian effects are present and that these effects must be
taken into account when attempting to model observable features of folding
proteins.

## Methods

In Section 1 of the SI document, we present details for the molecular dynamics
simulations, including various relevant simulation parameters. Additionally,
we include a range of measured time-scales and other extracted quantities.

The fraction of native contacts: For each protein, we project the back-bone Cα
atomic positions from the all-atom trajectories onto the fraction of native
contacts reaction coordinate $Q$, evaluated with a soft cut-off potential
[28]. The evaluation of $Q(t)$ requires a reference state, which we take to be
the native state for each trajectory. To evaluate the native state, we follow
previous implementations and select from amongst the member states of the
trajectories (as opposed to using, for example, the PDB entry). The approach
is similar to that used by Lindorff-Larsen et. al. [22] and Best et. al. [28],
which follows from [51]. Briefly, we sample a subset of evenly spaced states
from the full trajectory. For each pair of states, we calculate the
corresponding root-mean-squared deviation (RMSD) between the two states. If
the RMSD between two states is less than 0.2 nm, then we place the pair into a
list. We assign the state that has the most listed pairs satisfying the RMSD
condition as the native state for a given protein. Note that for proteins with
more than one independent trajectory segments we select a single native state
from amongst all trajectory segments, which we then use for all segments. In
the native state, we define all Cα pairs that are separated by at least 5
residues in the primary sequence and which are separated by less than 0.9 nm
in Cartesian distance, as the native contacts. Each protein will have
$N_{\rm{nc}}$ native contacts. $\textbf{s}_{ij}^{0}$ are the separation
vectors for all native contact pairs in the native state, which have
magnitudes $s_{ij}^{0}=\sqrt{\textbf{s}_{ij}^{0}\cdot\textbf{s}_{ij}^{0}}$.
$\textbf{s}_{ij}(t)$ are the separation vectors for all native contact pairs
at each time, with magnitudes
$s_{ij}(t)=\sqrt{\textbf{s}_{ij}(t)\cdot\textbf{s}_{ij}(t)}$. This gives the
fraction of native contacts we deem to be in contact at time t:

$Q(t)=\frac{1}{{N_{\rm{nc}}}}\sum^{N_{\rm{nc}}}_{i<j}\frac{1}{1+\text{e}^{\beta(s_{ij}(t)-\gamma
s_{ij}^{0})}}.$ (M1)

Here, we set the parameters such that $\beta=30$ nm-1 and $\gamma=1.6$.

Non-Markovian correction factor $\xi_{\text{noMar}}$: From Kappler et. al.
[13] and Lavacchi et. al. [14], we describe multi-exponential memory dependent
barrier crossing times as a sum of contributions from $M$ over-damped
contributions $\\{\tau^{i}_{\rm{OD}}\\}$, and $M$ energy-diffusion
contributions $\\{\tau^{i}_{\rm{ED}}\\}$, where $i=1,2,...,M$. The individual
contributions are defined as follows. For the over-damped contributions, we
have:

$\begin{split}\tau^{i}_{\rm{OD}}&=\tau_{\text{D}}\frac{\gamma_{i}}{\gamma}\frac{\text{e}^{\beta
U_{0}}}{\beta U_{0}}\\\
&\qquad{\times}\Bigg{[}\frac{\pi}{2\sqrt{2}}\frac{1}{1+10\beta
U_{0}\tau_{i}/\tau_{\text{D}}}+\sqrt{\beta
U_{0}\frac{\tau_{\rm{m}}}{\tau_{\text{D}}}}\Bigg{]},\end{split}$ (M2)

and for the energy-diffusion contributions:

$\begin{split}\tau^{i}_{\rm{ED}}&=\tau_{\text{D}}\frac{\gamma}{\gamma_{i}}\frac{\text{e}^{\beta
U_{0}}}{\beta U_{0}}\\\
&\qquad{\times}\Bigg{[}\frac{\tau_{\rm{m}}}{{\tau_{\text{D}}}}+4\beta
U_{0}\bigg{(}\frac{\tau_{i}}{\tau_{\text{D}}}\bigg{)}^{2}+\sqrt{\beta
U_{0}\frac{\tau_{\rm{m}}}{\tau_{\text{D}}}}\Bigg{]}.\end{split}$ (M3)

Here, $\beta=1/k_{\rm{B}}T$. We combine Eqs. M2 and M3 such that the predicted
mean passage times are given by:

$\begin{split}&\tau^{\text{H}}_{\rm{MFP}}(\tau_{\rm{D}},U_{0},\\{\gamma_{i}\\},\\{\tau_{i}\\},\tau_{\rm{m}})\\\
&\qquad\qquad\qquad\qquad=\sum_{i=1}^{M}\tau^{i}_{\rm{OD}}+\Bigg{[}\sum_{i=1}^{M}\frac{1}{\tau^{i}_{\rm{ED}}}\Bigg{]}^{-1}.\end{split}$
(M4)

$\\{\gamma_{i}\\}$ and $\\{\tau_{i}\\}$ are the sets of $M$ amplitudes and
time-scales that appear in Eqs. M2 and M3. The total friction $\gamma$ that
appears in Eqs. M2 and M3 is account for since
$\gamma=\sum_{i=1}^{M}\gamma_{i}$. The over-damped, Markovian limit is
achieved by letting all memory time-scales and inertial times equal to 0. The
limit is therefore given by:

$\tau_{\rm{MFP}}^{\text{H}}(\tau_{\rm{D}},U_{0},\\{\gamma_{i}\\},0,0)=\tau_{\text{D}}\pi
U_{0}\text{e}^{\beta U_{0}}/2\beta\sqrt{2},$ (M5)

leading to the non-Markovian correction factor $\xi_{\text{noMar}}$ and the
non-Markovian barrier crossing time $\tau_{\rm{MFP}}^{\rm{noMar}}$, as given
by Eq. 4 and Eq. 5, respectively.

## Acknowledgements

We are extremely grateful to the group of David E. Shaw for providing us with
the MD protein trajectories, which we have analysed in this paper. We are
grateful to the physics-department HPC services at the Freie University of
Berlin, for their generous provision of compute-time. This project was funded
by the European Research Council (ERC) Advanced Grant 835117 NoMaMemo.

## References

  * Bryngelson and Wolynes [1989] J. D. Bryngelson and P. G. Wolynes, The Journal of Physical Chemistry 93, 6902 (1989).
  * Bryngelson  _et al._ [1995] J. D. Bryngelson, J. N. Onuchic, N. D. Socci, and P. G. Wolynes, Proteins: Structure 21 (1995).
  * Dill and Chan [1997] K. A. Dill and H. S. Chan, Nature Structural Biology 4, 10 (1997).
  * Schuler and Eaton [2008] B. Schuler and W. A. Eaton, Current Opinion in Structural Biology 18, 16 (2008).
  * Hinczewski  _et al._ [2013] M. Hinczewski, J. C. M. Gebhardt, M. Rief, and D. Thirumalai, Proceedings of the National Academy of Sciences 110, 4500 (2013).
  * Plotkin and Wolynes [2003] S. S. Plotkin and P. G. Wolynes, Proceedings of the National Academy of Sciences 100, 4417 (2003).
  * Arrhenius [1889] S. Arrhenius, Zeitschrift für Physikalische Chemie 4U, 226 (1889).
  * Eyring [1935] H. Eyring, The Journal of Chemical Physics 3, 107 (1935).
  * Kramers [1940] H. A. Kramers, Physica 7, 284 (1940).
  * Mel’nikov and Meshkov [1986] V. I. Mel’nikov and S. V. Meshkov, The Journal of Chemical Physics 85, 1018 (1986).
  * Grote and Hynes [1980] R. F. Grote and J. T. Hynes, The Journal of Chemical Physics 73, 2715 (1980).
  * Kappler  _et al._ [2018] J. Kappler, J. O. Daldrop, F. N. Brünig, M. D. Boehle, and R. R. Netz, The Journal of Chemical Physics 148, 14903 (2018).
  * Kappler  _et al._ [2019] J. Kappler, V. B. Hinrichsen, and R. R. Netz, The European Physical Journal E 42, 119 (2019).
  * Lavacchi  _et al._ [2020] L. Lavacchi, J. Kappler, and R. R. Netz, EPL 131 (2020).
  * Echeverria  _et al._ [2014] I. Echeverria, D. E. Makarov, and G. A. Papoian, Journal of the American Chemical Society 136, 8708 (2014).
  * Schulz  _et al._ [2012] J. C. F. Schulz, L. Schmidt, R. B. Best, J. Dzubiella, and R. R. Netz, Journal of the American Chemical Society 134, 6273 (2012).
  * de Sancho  _et al._ [2014] D. de Sancho, A. Sirur, and R. B. Best, Nature Communications 5, 4307 (2014).
  * Lange and Grubmüller [2006] O. F. Lange and H. Grubmüller, The Journal of Chemical Physics 124, 214903 (2006).
  * Best and Hummer [2006] R. B. Best and G. Hummer, Physical Review Letters 96, 228104 (2006).
  * Best and Hummer [2010] R. B. Best and G. Hummer, Proceedings of the National Academy of Sciences 107, 1088 LP (2010).
  * Hinczewski  _et al._ [2010] M. Hinczewski, Y. von Hansen, J. Dzubiella, and R. R. Netz, The Journal of Chemical Physics 132, 245103 (2010).
  * Lindorff-Larsen  _et al._ [2011] K. Lindorff-Larsen, S. Piana, R. O. Dror, and D. E. Shaw, Science 334, 517 LP (2011).
  * Shaw  _et al._ [2009] D. E. Shaw, R. O. Dror, J. K. Salmon, J. P. Grossman, K. M. Mackenzie, J. A. Bank, C. Young, M. M. Deneroff, B. Batson, K. J. Bowers, E. Chow, M. P. Eastwood, D. J. Ierardi, J. L. Klepeis, J. S. Kuskin, R. H. Larson, K. Lindorff-Larsen, P. Maragakis, M. A. Moraes, S. Piana, Y. Shan, and B. Towles, in  _Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis_ (2009) pp. 1–11.
  * Shaw  _et al._ [2010] D. E. Shaw, P. Maragakis, K. Lindorff-Larsen, S. Piana, R. O. Dror, M. P. Eastwood, J. A. Bank, J. M. Jumper, J. K. Salmon, Y. Shan, and W. Wriggers, Science (New York, N.Y.) 330, 341 (2010).
  * Zwanzig [1961] R. Zwanzig, Physical Review 124, 983 (1961).
  * Mori [1965] H. Mori, Progress of Theoretical Physics 33, 423 (1965).
  * Nymeyer  _et al._ [1998] H. Nymeyer, A. E. García, and J. N. Onuchic, Proceedings of the National Academy of Sciences 95, 5921 (1998).
  * Best  _et al._ [2013] R. B. Best, G. Hummer, and W. A. Eaton, Proceedings of the National Academy of Sciences 110, 17874 LP (2013).
  * Ayaz  _et al._ [2022] C. Ayaz, L. Scalfi, B. A. Dalton, and R. R. Netz, Physical Review E 105, 54138 (2022).
  * Straub  _et al._ [1987] J. E. Straub, M. Borkovec, and B. J. Berne, The Journal of Physical Chemistry 91, 4995 (1987).
  * Daldrop  _et al._ [2018] J. O. Daldrop, J. Kappler, F. N. Brünig, and R. R. Netz, Proceedings of the National Academy of Sciences 115, 5169 (2018).
  * Darve  _et al._ [2009] E. Darve, J. Solomon, and A. Kia, Proceedings of the National Academy of Sciences 106, 10884 (2009).
  * Ayaz  _et al._ [2021] C. Ayaz, L. Tepper, F. N. Brünig, J. Kappler, J. O. Daldrop, and R. R. Netz, Proceedings of the National Academy of Sciences 118, e2023856118 (2021).
  * Kowalik  _et al._ [2019] B. Kowalik, J. O. Daldrop, J. Kappler, J. C. F. Schulz, A. Schlaich, and R. R. Netz, Physical Review E 100, 12126 (2019).
  * Einert  _et al._ [2011] T. R. Einert, C. E. Sing, A. Alexander-Katz, and R. R. Netz, The European Physical Journal. E, Soft Matter 34, 1 (2011).
  * de Gennes [1971] P. G. de Gennes, The Journal of Chemical Physics 55, 572 (1971).
  * Chung and Eaton [2013] H. S. Chung and W. A. Eaton, Nature 502, 685 (2013).
  * Sung  _et al._ [2012] C. H. Sung, M. Kevin, L. J. M., and E. W. A., Science 335, 981 (2012).
  * Chung  _et al._ [2015] H. S. Chung, S. Piana-Agostinetti, D. E. Shaw, and W. A. Eaton, Science (New York, N.Y.) 349, 1504 (2015).
  * Krishna  _et al._ [2016] N. Krishna, F. D. A. N., D. D. R., Y. Hao, W. Feng, and W. M. T., Science 352, 239 (2016).
  * Satija  _et al._ [2017] R. Satija, A. Das, and D. E. Makarov, The Journal of Chemical Physics 147, 152707 (2017).
  * Daldrop  _et al._ [2016] J. O. Daldrop, W. K. Kim, and R. R. Netz, EPL (Europhysics Letters) 113, 18004 (2016).
  * Hummer [2003] G. Hummer, The Journal of Chemical Physics 120, 516 (2003).
  * Laleman  _et al._ [2017] M. Laleman, E. Carlon, and H. Orland, The Journal of Chemical Physics 147, 214103 (2017).
  * Carlon  _et al._ [2018] E. Carlon, H. Orland, T. Sakaue, and C. Vanderzande, The Journal of Physical Chemistry B 122, 11186 (2018).
  * Sali  _et al._ [1994] A. Sali, E. Shakhnovich, and M. Karplus, Nature 369, 248 (1994).
  * Du  _et al._ [1998] R. Du, V. S. Pande, A. Y. Grosberg, T. Tanaka, and E. S. Shakhnovich, The Journal of Chemical Physics 108, 334 (1998).
  * Best and Hummer [2005] R. B. Best and G. Hummer, Proceedings of the National Academy of Sciences 102, 6732 LP (2005).
  * Pirchi  _et al._ [2016] M. Pirchi, R. Tsukanov, R. Khamis, T. E. Tomov, Y. Berger, D. C. Khara, H. Volkov, G. Haran, and E. Nir, The Journal of Physical Chemistry B 120, 13065 (2016).
  * Aviram  _et al._ [2018] H. Y. Aviram, M. Pirchi, H. Mazal, Y. Barak, I. Riven, and G. Haran, Proceedings of the National Academy of Sciences 115, 3243 (2018).
  * Daura  _et al._ [1999] X. Daura, K. Gademann, B. Jaun, D. Seebach, W. F. van Gunsteren, and A. E. Mark, Angewandte Chemie International Edition 38, 236 (1999).