# Impact of hydraulic tortuosity on micro/nanoporous flow

Shiwani Singh<EMAIL_ADDRESS>Mathematics Institute, University of
Warwick, Coventry CV4 7AL, United Kingdom

###### Abstract

Using the porous structures made up of homogeneously arranged solid obstacles,
we examine the effects of rarefaction on the hydraulic tortuosity in the slip
and early transition flow regimes via extended lattice Boltzmann method. We
observed that modification in either the obstacle’s arrangement or the
porosity led to a power-law relation between the porosity-tortuosity. Along
with this, we also found that in the slip flow regime, the exponent of this
relation contains the effect of finite Knudsen number (Kn). In addition, we
observed that on properly scaling Kn with porosity and hydraulic tortuosity, a
generalized correlation can be obtained for apparent permeability.

##  I Introduction

Numerous scientific and engineering applications, including water percolating
through soils, gas transportation, CO2 sequestration, oil extraction with or
without polymer flooding and a host of others, can benefit from an
understanding of the transportation mechanisms and fluid flow across porous
media bear1972dynamics ; ghanbarian2013tortuosity ; metz2005ipcc ;
lake1989enhanced . Porous media are made up of a solid matrix of material that
is filled by a network of void spaces (pores) containing fluid and are
connected by throats that are significantly smaller in size. One of the main
physical properties of interest is the permeability of a porous material which
measures its ability to allow fluids (gas or liquid) to flow through it. Apart
from porosity which is defined as the percentage of a porous sample that is
occupied by pore space, permeability also depends on the geometry and
structure of the pores carman1937fluid ; bear1972dynamics ; dullien2012porous
.

Depending on the geometry and the location of the pores, the actual path taken
by the fluid can be very complicated or tortuous. Therefore, a parameter,
hydraulic tortuosity (T), was introduced to take care of complicated transport
paths in a comprehensive manner bear1972dynamics ; koponen1996tortuous .
Hydraulic tortuosity can be understood as the ratio of the average length of
true flow routes to the system’s length in the direction of the macroscopic
flow. The optimal way to calculate this ratio would be to take the mean flow
length from a weighted average of the streamline, however if the geometry is
too complicated, it would not be feasible matyka2008tortuosity ;
ghanbarian2013tortuosity . To overcome this limitation, Koponen et al.
koponen1996tortuous  and Duda et al. duda2011hydraulic  proposed a
straightforward formula for T as a ratio between the mean of total fluid
velocity and the mean of component of the velocity along the external force
direction, which is extremely helpful in the situation of arbitrary geometry.

Furthermore, the emergence of unconventional energy sources, like ultra-tight
gas reservoirs within shale rocks, have shown great potential towards
mitigating the world energy crisis javadpour2007nanoscale ;
freeman2011numerical . Shale rocks are highly tortuous and are made up of
fine-grained material which contains pores in the nanoscale size range. At
this scale where the mean free path of gas molecules becomes equal to or
greater than the characteristic flow length within ultra-tight rocks, the
rarefaction effects starts to emerge. The Knudsen number Kn, or the ratio of
the mean free route of gas molecules to the typical flow length, indicates the
degree of rarefaction. The fluid behavior can be separated into four primary
categories based on the Kn value: continuum flow regime with Kn $<$0.001, slip
flow regime with $0.001<{\rm Kn}<0.1$, transition flow regime with $0.1<{\rm
Kn}<10$, and free molecular flow regime with Kn $>10$. Contrary to what is
predicted by Darcy’s law in continuum flow regime, rarefaction effect cause
gas permeability (apparent) to increase as the pore size decreases. Kinkenberg
was the first to claim that this increase is caused due to the rise in gas
slippage at the solid-fluid interface klinkenberg1941permeability . Therefore,
to precisely predict the reservoir’s production capacity and longevity, it
becomes unquestionably crucial to investigate the effect of rarefaction on
physical parameters like permeability and tortuosity of ultra-tight porous
media.

The goal of this research is to provide a deeper understanding of the gas
transport characteristics and Kn dependency of various physical properties of
porous media, including tortuosity, in carefully designed porous media. The
set-up is designed in a simple way where the circular obstacles are arranged
homogeneously between two parallel plates. The porosity is varied by changing
the diameter of the obstacle and the tortuosity is altered by changing the
location of next nearest obstacle. As a simulation technique, we used the
lattice Boltzmann (LB) method which has not only proven to be a useful tool
for simulating Newtonian continuum hydrodynamics benzi1992lattice ;
chen1998lattice ; succi2001lattice ; aidun2010lattice ; kruger2017lattice  but
also been successfully extended for flows beyond Navier-Stokes in the past few
years particularly to non-equilibrium ( finite Kn flows). The extended LB
method utilizes either a regularization procedure zhou2006simulation ;
zhang2006efficient ; latt2006lattice  or an appropriate multi-relaxation time
model tao2015boundary ; guo2014generalized  in combination with kinetic
boundary condition ansumali2002kinetic ; montessori2015lattice ;
singh2017impact ; singh2017influence . We also conducted a thorough parametric
study with the goal of determining a consistent way to account for the impact
of porosity and tortuosity on gas permeability.

Starting from a brief overview of tortuosity and its evaluation technique in
Section II, the rest of the paper is organized as follow: In Section III, the
representative porous media set-up is detailed and a reference has been made
to simulation technique, the extended lattice Boltzmann method, which is
detailed later in the Appendix A. This followed by the Kn depended flow
investigation in Section IV where starting by investigating the local flow
profile at various physical parameter like porosity and porous arrangement, we
studied the effect of Kn on tortuosity-porosity relation in Section IV.1.
Further, the effect of tortuosity and porosity on gas permeability was studied
in Section IV.2. Finally, the work is summarized and some future aspects of
the work are discussed in Section V.

##  II  Tortuosity

Hydraulic tortuosity (T) is the ratio of elongation of fluid streamlines due
to the presence of obstacle (porous media) to the system size in the case of
free flow matyka2008tortuosity ; ghanbarian2013tortuosity . Therefore, if
$\lambda$ is the mean distance covered by fluid element and L is the system
size in the direction of flow, the hydraulic tortuosity is defined as

$T=\frac{\lambda}{L}.$ (1)

By this definition, hydraulic tortuosity is always greater than or equal to 1
$(T\geq 1)$. It means in a plane channel flow, $T=1$ since streamlines face no
hindrance in the absence of of porous material. The value of $T$ rises when
tortuosity grows because fluid has to travel farther through porous media.
Since T is defined as deviation of fluid path, it can be calculated using
velocity field. Some of the methods used to calculate $T$ focuses on the
calculation of weighted averages of discrete streamline knackstedt1994direct ;
zhang1995direct . Most of these methods will be difficult to apply in complex
geometries. However, in the works of Koponen et al. koponen1996tortuous  and
Duda et al. duda2011hydraulic , the authors came up with a simple approach for
the incompressible and non-reentrant flow to calculates $T$ directly from the
velocity field by averaging its components in the following manner:

$T=\frac{\langle u\rangle}{\langle u_{x}\rangle}$ (2)

where $\langle u\rangle$ is the average of the magnitude of fluid velocity and
$\langle u_{x}\rangle$ is the average of its component along the direction of
flow. Hereafter, in the present study, we will be using the above mentioned
definition (Eq. (2)) which is now being routinely used to define $T$
icardi2014pore ; matyka2016power ; muljadi2016impact ; el2019pore ;
graczyk2020predicting . A detailed pedagogical review about the approach used
in the Ref.duda2011hydraulic  can be found in the Ref. matyka2012calculate .

##  III  Physical porous media

In the present study, we used simple two-dimensional homogeneous geometries
constructed by placing circular obstacle in seven rows and seven column
between two parallel plates to study the effect of tortuosity (see Fig. 1).
The lattice Boltzmann (LB) method is used to solve the flow equations. In the
recent past, the LB method has emerged as an powerful tool to solve the
continuum hydrodynamics benzi1992lattice ; succi2001lattice  and with some
physics inspired refinements has shown to capture fluids dynamics beyond
Navier-Stokes equationsmontessori2015lattice ; singh2017impact ;
mohammed2021lattice . The details of this refined method is presented in the
Appendix A.

In the aforementioned geometrical set-up, the lattice node which falls on or
inside solid is marked as ‘$\rm node_{S}$’ and the ones occupied by fluid are
marked $`\rm node_{F}$’. Hence, the calculation of the porosity( $\phi$),
defined as the volume (area in two-dimension) of void to the total volume
(area), is straight forward and is given as

$\phi=\frac{\rm node_{F}}{\rm node_{F}+\rm node_{S}}.$ (3)

(a) $\theta_{0}$, $\phi=0.75$

(b) $\theta_{6}$, $\phi=0.75$

(c) $\theta_{6}$, $\phi=0.90$

Figure 1: Schematic of representative porous media designed as an array of
circular obstacle in two-dimension

At this point, it is worth introducing ${\rm Kn}$ which is the ratio of the
molecular mean free path with respect to character macroscopic length. The Kn
can also be represented in terms of kinematic viscosity, $\eta$, as ${\rm
Kn}=\eta/(dc_{s})$ where we chose $d$ to be the smallest pore-throat diameter
present as red double arrow in the Fig. 1, and $c_{s}$ is the speed of sound
which is explained later in the Appendix A where LB method is elaborated
alongside with regularization mechanism which filters out the non-hydrodynamic
moments and the kinetic boundary condition (KBC) which is based on diffusively
reflecting wall. In the recent past, it has been observed that regularization
and KBC are the crucial ingredients required in the standard LB for the
simulation of finite Kn flow montessori2015lattice ; singh2017impact . The
simulations were performed with 500 grid points in each direction which
resulted in 24-38 number grid points representing each pore throat.

Before studying the Kn dependent rarefaction effects, we focused in the
continuum regime to establish the relation between the porosity and
tortuosity. Consequently, starting from uniform geometry where distance
between the next nearest circle placed in row is the same as that of the one
placed in column (Fig. 1(a)), the simplest way to alter tortuosity is to
change the arrangement by placing the next nearest circle at a distance
defined by the angle $\theta$ as shown in Fig. 1(b). This allowed us to do a
controlled study and have insight of how tortuosity increases/decreases with
the arrangement of obstacles in the media.

Figure 2: Tortuosity as a function of porosity for various arrangement in
porous media. Solid lines are the best fit to $T-1=p_{n}(1-\phi)^{\gamma_{n}}$
where $n(0,1,...,6)$ defines the alignment.

We choose the following values of alignments:
$[\theta_{0},\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5},\theta_{6}]=[0,5.5^{\circ},9.162^{\circ},10.45^{\circ},14.47^{\circ},17.88^{\circ},21.96^{\circ},25.82^{\circ}]$.
The streamlines passing through the pore throat, present nearly at the middle
of the domain, marked with a red double arrow shows that keeping the porosity
same, tortuosity appears to increase with $\theta$ (see Fig. 1(a) and Fig.
1(b)). In Fig. 1(c), we kept the alignment ($\theta$) same as that of Fig.
1(b) but increased the porosity which resulted in an obvious decrease in
tortuosity. Recall, as defined earlier, larger the average length of
streamlines, larger will be the tortuosity. To calculate the absolute
tortuosity, $T_{\infty}$, for each porous configuration in the limit of Kn
$\to$ 0, we extrapolated the tortuosity (T) calculated within the range from
${\rm Kn}=10^{-1}$ to ${\rm Kn}=10^{-3}$ and chose the value at ${\rm
Kn}=10^{-7}$ as $T_{\infty}$.

Fig. 2 represents the the absolute tortuosity ($T_{\infty}$) as function of
porosity for large to medium porosity ranging from 0.90 to 0.75 for all the
seven alignment ($\theta_{0}\textrm{ to }\theta_{6}$). Firstly, the figure
shows the obvious trend that with increase in porosity, tortuosity decrease
for every configuration. Secondly, the log-log plot between $T_{\infty}-1$ and
$1-\phi$ clearly shows a power-law behaviour with the exponent as
$[\gamma^{0},\gamma^{1},\gamma^{2},\gamma^{3},\gamma^{4},\gamma^{5},\gamma^{6}]=\left[0.1286,0.3153,0.4199,0.4902,0.5706,0.6530,0.7279\right]$
and constants as
$[p_{0},p_{1},p_{2},p_{3},p_{4},p_{5},p_{6}]=\left[0.0233,0.0412,0.0702,0.1084,0.1606,0.2243,0.2866\right]$
(see Fig. 2). This is an empirical relationship, however, is in agreement with
the two-dimensional flow with randomly distributed square obstacle as used by
Duda et. al duda2011hydraulic  for large porosity ($<0.8$) which also showed a
power-law behaviour as $T-1\sim(1-\phi)^{\gamma}$, however their exponent
($\gamma$) was $1/2$.

##  IV  Flow at finite Knudsen Number

Figure 3: Steady state velocity streamlines for porosity, $\phi=0.75$ and
$\phi=0.90$, porous arrangement , $\theta_{0}$ and $\theta_{6}$, and at
Kn=0.01, Kn=0.1 and Kn=1.0.

Tortuosity and permeability are two important properties that influence the
passage of fluid through a porous media. In the subsequent sections, we will
investigate the effect of rarefaction on these two properties. However,
firstly, it would also be interesting to observe the behavior of the local
velocity in different porous arrangement. Therefore, in Fig. 3, we plotted
steady state streamlines at different Kn for porosities, $\phi=0.75$ and
$\phi=0.90$, and at pore alignments, $\theta_{0}$ and $\theta_{6}$. The scale
for the magnitude of streamlines in all the cases has been kept the same to
show a relative difference between all the considered arrangements. From the
Fig. 3, we can make the following observations:

  * •

For a given porosity ($\phi$), the flow take more tortuous path as the
uniformity of the obstacles decreases (i.e. from $\theta_{0}$ to $\theta_{6}$)
at every Kn.

  * •

In the current configuration, when porosity rises, the pore throat widens,
facilitating fluid flow across the porous medium. As a result, the velocity
inside the pore-throat increased as the $\phi$ value increased.

  * •

However, due to Knudsen diffusion, there is a non-zero fluid velocity at solid
barriers, which causes the velocity inside the pore throat to grow with
increasing Kn in all situations.

###  IV.1  Effect of finite Kn on tortuosity

(a) $\theta_{1}$

(b) $\theta_{3}$

(c) $\theta_{6}$

Figure 4: Tortuosity as a function porosity at pore alignment $\theta_{1}$,
$\theta_{3}$ and $\theta_{6}$ with varying Kn. Dashed lines are the best fit
to the function $T-1=p_{n}(1-\phi)^{(\gamma_{n}+{\rm f}_{n}({\rm Kn}))}$ where
the value of $n$(0,1,2,..6) depends on the porous arrangement.

Hydraulic tortuosity (T) is the measure of average fluid streamline length in
a porous medium versus system length of obstacle-free flow, and as defined in
Eq. 1, can be calculated using appropriate averages of velocity field
component. This makes T a flow dependent observable. Therefore, any effects of
rarefaction that appears on velocity field such as, increment of slip velocity
with Kn, also have significant effect on T.

To investigate this effect, we computed the value of T as a function of
porosity, $\phi$, with different Kn for all the arrangement of porous medium,
$\theta_{0}$-$\theta_{6}$ which again shows a power law behavior, however with
a different exponent (Kn-dependent) as shown in Fig. 4. When Kn rises, the
slip velocity makes it easier for fluid to move through the pores, which
causes a concomitant decrease in the hydraulic tortuosity. This behavior is
clearly reflected in Fig. 4. In order to fully contain the Kn-dependency of
hydraulic tortuosity in the exponent of the power-law behavior, we plotted the
$T-1$ and $1-\phi$ with the best fitted lines to the function
$T-1=p_{n}(1-\phi)^{(\gamma_{n}+{\rm f}_{n}({\rm Kn}))}$ in Fig. 4 where $n$
(0,1,2,..6) dictates the pore arrangement. Here, the values of $\gamma_{n}$
and $p_{n}$ are the same to those assessed in Section III (Fig. 2). Changing
Kn while maintaining the same configuration demonstrates that ${\rm f}_{n}$ is
an increasing function of Kn because as Kn rises, the exponent of the power
law rises as well. Fig. 4 also shows that the exponent in the early-transition
phase (bottom few lines in each figure) appears to change less profoundly with
changing Kn as the non-uniformity between the pores increases from Fig. 4(a)
($\theta_{1}$) to Fig. 4(b) ($\theta_{3}$) and eventually to Fig. 4(c)
($\theta_{6}$). This shows that when Kn is increased at high non-uniformity,
there is less of an effect on physical attributes. In other words, the effect
of rarefaction in early transition regime is diminished by higher non-
uniformity. But this can also be a result of a technical problem associated
with the method. It has been shown in earlier research that a lattice
Boltzmann model with regularization and a diffuse wall kinetic boundary
condition (as used in the current study) performs admirably well in the slip
flow regime and near to the early-transition flow flow. However, to take into
account the flow beyond this limit, one must utilize a higher order lattice,
according to Refmontessori2015lattice . This aspect of scheme needs further
investigation, thus it is left for future study.

Figure 5: The Kn dependent addition in the exponent of power-law behavior,
$f_{n}$, for various Kn. Here $n$ (0,1,…,6) describes the alignment of porous
arrangement. The inset magnifies the slip flow regime which shows that $f({\rm
Kn})\sim{\rm Kn}^{0.6}$.

Furthermore, to investigate the behavior of ${\rm f}_{n}$ with respect to Kn,
we plotted ${\rm f}_{n}$ as a function of Kn (see Fig. 5) which firstly
indicate that ${\rm f}_{n}$ is indeed an increasing function of Kn for all the
porous arrangements. However, in the slip regime ($0.001<{\rm Kn}<0.1$), ${\rm
f}_{n}\sim{\rm Kn}^{0.6}$ (see the inset of Fig. 5). This suggests that there
exists a consistent dependence of Kn on the tortuosity-porosity relation at
least in the slip flow regime for the porous setup used here as,
$T-1\sim(1-\phi)^{(\gamma_{n}+{\rm Kn}^{0.6})}$.

###  IV.2  Effect of finite Kn on permeability

One of the physical quantity of interest for the flow inside highly permeable
porous media is the absolute permeability ($\kappa_{\infty}$) which is
obtained by calculating the ﬂux ($Q$) at different pressure drops (or by
adjusting the body force values) in the following manner:

$\kappa_{\infty}=\frac{Q\eta}{\rho g}$ (4)

where $\eta$ is the dynamic viscosity, g is the body force and $\rho$ is the
density of the fluid. Permeability is a crucial element in determining the
transport capacity of porous media. However, in unconventional reservoirs such
a ultra-tight pores of shale rock, despite the presumption of absolute
permeability being very low, experiments observed that apparent gas
permeability (AGP) ($\kappa$), is much higher than $\kappa_{\infty}$.
Furthermore, the idea that the rise in gas slippage at the solid-fluid
interface is to blame for the increased permeability was first put forward by
Kinkenberg who suggested the permeability correction factor (PCF), which is
defined as the ratio of apparent permeability ($\kappa$) to absolute
permeability ($\kappa_{\infty}$), to be a linear function of Kn as:

$\frac{\kappa}{\kappa_{\infty}}=1+4{\rm Kn}.$ (5)

Beskok and Karniadakis beskok1999report  further proposed a second-order
correlation that can be used to describe all four fluid flow regimes and is
given as:

$\frac{\kappa}{\kappa_{\infty}}=[1+\alpha({\rm Kn}){\rm
Kn}]\left(1+\frac{4{\rm Kn}}{1-b{\rm Kn}}\right)$ (6)

where slip coefﬁcient $b$ equals -1 for slip ﬂow, and $\alpha({\rm Kn})$ is
the rarefaction coefﬁcient. The expression for $\alpha({\rm Kn})$ is somewhat
complex, but Civan civan2010effective  later suggested one (Beskok &
Karniadakis-Civan’s correlation) that is considerably more straightforward:

$\alpha({\rm Kn})=\frac{1.358}{1+0.170{\rm Kn}^{-0.4348}}.$ (7)

In addition, Civan civan2010effective  suggested that in the slip flow regime,
$\alpha({\rm Kn})$ can also be neglected making PCF take the following form:

$\frac{\kappa}{\kappa_{\infty}}=\left(1+\frac{4{\rm Kn}}{1+{\rm Kn}}\right).$
(8)

Since up-scaled equations that predict the gas output and longevity of gas
wells heavily rely on the AGP at representative elementary volume (REV) scale,
it is essential to accurately evaluate AGP. For this reason, we first
presented the permeability correction factor (PCF) of different porous medium
by varying the alignment at fixed porosity (Fig. 6) and then varied the
porosity while keeping the alignment unchanged (Fig. 7). To calculate the
absolute permeability ($\kappa_{\infty}$) for each set-up, we extrapolated the
permeability calculated within the range from Kn = $10^{-3}$ to Kn = $10^{-1}$
and chose the value at Kn=$10^{-7}$ as $\kappa_{\infty}$.

(a) $\phi=0.75$

(b) $\phi=0.83$

(c) $\phi=0.90$

Figure 6:  The permeability correction factor (PCF) as a function of Kn at
fixed porosity and varying alignments.

When the porosity is low ($\phi=0.75$) as shown in Fig. 6(a), the setup with
more non-uniformly distributed porous media shows a slightly higher value of
PCF over the range of Kn. However, as the porosity is increased to higher
values as one shown in Fig. 6(b) ($\phi=0.83$) and Fig. 6(c)($\phi=0.90$), the
PCF starts to overlap on to each other for all the alignments.

(a) $\theta_{0}$

(b) $\theta_{3}$

(c) $\theta_{6}$

Figure 7:  The permeability correction factor (PCF) as a function of Kn at
fixed alignment and varying porosity.

Likewise, in Fig. 7, we kept the alignment (uniformity/non-uniformity) of
distributed porous media fixed and varied the porosity. In this case, we
observed that as Kn is increased (especially in the early transition regime)
in all the alignments, the PCF for low porosity has a higher value in
comparison to the one with higher porosity. Additionally, this difference
becomes more prominent as the dispersion of porous media is further increased
from $\theta_{3}$ to $\theta_{6}$. One reason for this behavior at the
macroscale level could be that in less porous media the disturbance caused by
one barrier propagates over a smaller distance and interacts with other
obstacles before getting completely dispersed. The hydrodynamic disturbance
may, however, tend to fade out or average out in highly porous media before
encountering another barrier. This causes different kinds of fluid-solid
interactions in addition to slip effect, which could change the magnitude of
the velocity field. At the molecular level, large pores with high porosity
have a lower rarefaction impact than smaller pores with low porosity because
after colliding inside one pore’s throat, the gas molecule must travel a
greater distance before colliding with another obstruction in highly porous
medium as compared one with small porosity. Overall, figures 6 and 7 led us to
the conclusion that, in the current porous setting, porosity has a greater
impact on PCF than does non-uniformity (alignment) of porous material, which
is one of the factors in determining the tortuosity of the material.

Figure 8: The permeability correction factor (PCF) as a function of an
effective ${\rm Kn}^{\star}$ at randomly selected porous media configuration.
The first-order correlation by Klinkenberg (Eq. (5))
klinkenberg1941permeability , the second-order correlation given by Beskok &
Karniadakis-Civan (Eq. (6))beskok1999report ; civan2010effective  and the
correlation in the slip flow regime suggested by Civan (Eq.
(8))civan2010effective  are also presented as a reference.

In recent studies, the direct simulation Monte Carlo (DSMC) approach was
employed to demonstrate that even for complicated porous structures, the
apparent permeability of a porous media can be predicted through the
Klinkenberg correlations using fundamental and measurable descriptors of the
pore structure wu2016non ; wu2017apparent ; yang2018investigation ; su2020gsis
. This can be accomplished by scaling Kn with an appropriate function of
porosity and tortuosity, as suggested by Wu et al. wu2017apparent , or by
using an effective pore size that is determined by porosity and tortuosity
yang2018investigation , as it was noticed that using such arguments, the PCF
for different porosity and pore structures is found to collapse onto one
curve. In a similar manner, we identified an empirical function,
$g=\phi/T^{2}$, to determine an effective Knudsen number ($\rm Kn^{\star}$)
for the current collection of porous media, given as ${\rm Kn}^{\star}={\rm
Kn}/{g}$. This can be inferred indirectly as rescaling of the pore size using
parameters related to pore structure, porosity, and pore scale flow, hydraulic
tortuosity (T). In Fig. 8, the PCF is plotted with respect to the effective
Knudsen number ${\rm Kn}^{\star}$. The figure shows that for all the
geometries (chosen at random), the PCF lie almost entirely along a single
line. Figure 8 further demonstrates that the PCF and the three proposed
correlations- Klinkenberg (Eq. (5)) klinkenberg1941permeability , Beskok &
Karniadakis-Civan (Eq. (6))beskok1999report ; civan2010effective , and Civan
(Eq. (8))civan2010effective -are in fair agreement up until the slip flow
regime. However, in the early transition regime, the best agreement is with
Civan civan2010effective , which was initially suggested for the flow close to
slip regime. Figure 3 illustrates that a generalized correlation can be
obtained up until the early transition regime with an appropriate effective
Knudsen number.

##  V Conclusion and Outlook

We investigated the gas transport characteristics flowing across porous medium
containing micro- and nanopores for moderate to high porosities of
$0.75<\phi<0.90$ from the continuum to early-transition regimes in a
systematic manner using the extended lattice Boltzmann method with a focus on
two physical parameters of the flow namely, tortuosity and permeability. A
priori, it is not obvious that tortuosity and porosity in general have a
universal relationship. However, such correlation can arise, at least for some
porous material types. With this aim, we explored the gas transport in
homogeneously arranged porous media which showed an empirical power-law
behaviour between the two as $T-1\sim(1-\phi)^{\gamma}$. An further
investigation across various Kn shows that, in the slip flow regime, the Kn-
dependency alters the exponent of the power as
$T-1\sim(1-\phi)^{\gamma+f_{n}}$ with $f^{n}\sim{\rm Kn^{0.6}}$, giving a
generalized relation between the two.

In addition, we discovered that permeability correction factor (PCF) with
respect to Kn, which appears as a result of rarefaction, fall nearly on one
single line up until early transition regime, indicating a generalized
correlation, with appropriate scaling of Kn with parameters like porosity and
hydraulic tortuosity. Moreover, how the Kn should be properly scaled using
$\phi$ and $T$ in a heterogeneous and tortuous environment may also be crucial
to examine other phenomena such as the dispersive transport of a scalar
through complex media with pores size ranging from micro to nano-meters
meigel2022dispersive . It should be emphasized that the implementation of the
diffuse wall boundary condition, which requires wall-normals, was simple due
to the simplicity of the current 2D set-up formed by arranging circular
obstacles in a homogeneous pattern at different orientations. The distribution
of pores in a real rock, however, is highly heterogeneous and tortuous. In
future, we will explore the fluid flow in a realistic low-porosity geometry of
natural rock imaged by a multislice micro-CT scannerjiang2014changes ;
jiang2023upscaling  by utilizing a boundary condition that was more practical
for such a scenario and does not require calculating wall normals as proposed
in Ref. singh2017impact .

##  VI Acknowledgment

S.S. acknowledges the financial support by the Leverhulme Early Career
Fellowship (ECF-2019-100). S.S. also acknowledges the use of the Scientific
Computing Research Technology Platform and associated support services at the
University of Warwick, in the completion of this work.

##  Appendix A Lattice Boltzmann Method

The conventional LBM framework is based on the Boltzmann equation with single
relaxation approximation also known as (BGK) approximation bhatnagar1954model
. The discrete form of such a equation requires a set of discrete populations
$f=\\{f_{i}\\}$ and corresponding to which there exists a pre-defined discrete
velocities ${\boldsymbol{c}}_{i}$ $(i=1,\cdots,N)$ succi2001lattice ;
benzi1992lattice  and has the following form

$f_{i}({\boldsymbol{x}}+{\boldsymbol{c}}\Delta t,t+\Delta
t)=f_{i}({\boldsymbol{x}},t)+\Omega_{i}(f)+\Delta t{F}_{i}.$ (9)

Here, $F_{i}$ corresponds to $i$-th component of external force and the BGK
collision approximation, $\Omega_{i}(f)$, given as

$\Omega_{i}(f)=\frac{\Delta t}{\tau}(f_{i}^{\rm
eq}-f_{i}({\boldsymbol{x}},t)],$ (10)

dictates the relaxation of distribution function to an equilibrium Maxwell-
Boltzmann function, $f^{\rm eq}$, as given in Eq. (11), at the rate of
$\tau^{-1}$.

$\displaystyle\begin{split}{f}_{i}^{\rm
eq}&=w_{i}\rho\Biggl{[}1+\frac{{\boldsymbol{c}}_{i}\cdot{\boldsymbol{u}}}{\,c_{s}^{2}}+\frac{({\boldsymbol{c}}_{i}\cdot{\boldsymbol{u}})^{2}}{2\,\,c_{s}^{4}}-\frac{({\boldsymbol{u}}\cdot{\boldsymbol{u}})}{2\,\,c_{s}^{2}}\Biggr{]}.\end{split}$
(11)

The two-dimensional model (D2Q9) which is chosen for the present study has the
following nine discrete velocities

${\boldsymbol{c}}_{i}=\begin{cases}(0,0)&\mbox{if }i=0\\\
\left(\cos{\frac{(i-1)\pi}{4}},\sin{\frac{(i-1)\pi}{4}}\right)&\mbox{if
}i=1,2,3,4\\\
\sqrt{2}\left(\cos{\frac{(i-1)\pi}{4}},\sin{\frac{(i-1)\pi}{4}}\right)&\mbox{if
}i=5,6,7,8,\end{cases}$ (12)

with the corresponding weights as

$w_{i}=\begin{cases}\frac{4}{9}&\mbox{for }i=0\\\ \frac{1}{9}&\mbox{for
}i=1,2,3,4\\\ \frac{1}{36}&\mbox{for }i=5,6,7,8.\end{cases}$ (13)

The lattice sound speed $c_{s}$ that appeared in the Eq. (11) is related to
the magnitude of discrete velocity as $c^{2}=3c_{s}^{2}$. The relevant
hydrodynamic macroscopic moments, like density ($\rho$), momentum density
$(\rho{\boldsymbol{u}})$ and momentum flux (${\boldsymbol{\Pi}}$) can be
obtained by linear weighted sums as
$\rho({\boldsymbol{x}},t)=\sum_{i}f_{i},\quad\rho({\boldsymbol{x}},t){\boldsymbol{u}}({\boldsymbol{x}},t)=\sum_{i}f_{i}{\boldsymbol{c}}_{i},\quad\textrm{and}\quad{\boldsymbol{\Pi}}({\boldsymbol{x}},t)=\sum_{i}f_{i}({\boldsymbol{c}}_{i}{\boldsymbol{c}}_{i}-c_{s}^{2}{\boldsymbol{\delta}}),$
respectively with ${\boldsymbol{\delta}}$ being the identity matrix. Finally,
the term $F_{i}$ in Eq. (9), which is the $i^{\rm th}$ component of the body
force, is given as guo2002discrete :

$F_{i}=w_{i}\rho\left[\frac{{\boldsymbol{g}}\cdot{\boldsymbol{c}}_{i}}{c_{s}^{2}}+\frac{({\boldsymbol{g}}{\boldsymbol{u}}+{\boldsymbol{u}}{\boldsymbol{g}})}{2c_{s}^{2}}{\boldsymbol{:}}({\boldsymbol{c}}_{i}{\boldsymbol{c}}_{i}-c_{s}^{2}{\boldsymbol{\delta}})\right],$
(14)

where ${\boldsymbol{g}}$ is the constant acceleration vector. In the
subsequent section, we will briefly discuss a regularization scheme which was
introduced to filter out the nonphysical effect from finite Kn flow.

###  A.1 Regularization Scheme

Initially proposed to resolve the issue of stability of high viscous flows,
the regularization process, as introduced by Chen and co-workers
zhou2006simulation ; zhang2006efficient  and Latt and Chopard latt2006lattice
, turns out be one of the major ingredient for the finite Kn flow
montessori2015lattice  in the lattice Boltzmann framework. In regularization
process, non-hydrodynamic (ghost) modes are filtered out by dividing the post
streaming distribution function into two parts as:

$f_{i}=f_{i}^{\rm eq}+f_{i}^{\rm neq}.$ (15)

The information about the hydrodynamic modes is contained in $f_{i}^{\rm eq}$
and to remove the information of non-hydrodynamic modes from $f_{i}^{\rm
neq}$, it is converted into a new distribution function $f^{\rm Reg}_{i}$
which is then defined in terms of hydrodynamic moments
$(\rho,{\boldsymbol{u}},{\boldsymbol{\Pi}})$ and has the following discrete
form:

$f^{\rm
Reg}_{i}=\frac{w_{i}}{2c_{s}^{4}}({\boldsymbol{c}}_{i}{\boldsymbol{c}}_{i}-c_{s}^{2}{\boldsymbol{\delta}}):{\boldsymbol{\Pi}}^{\rm
neq},$ (16)

where ${\boldsymbol{\Pi}}^{\rm neq}$ is the non-equilibrium part of momentum
flux and is given as ${\boldsymbol{\Pi}}^{\rm neq}=\sum_{i}f^{\rm
neq}({\boldsymbol{c}}_{i}{\boldsymbol{c}}_{i}-c_{s}^{2}{\boldsymbol{\delta}})$.
The streaming step of lattice Boltzmann method takes the following form after
the regularization process:

$f_{i}({\boldsymbol{x}}+{\boldsymbol{c}}\Delta t,t+\Delta t)=f_{i}^{\rm
eq}+\left(1-\frac{\Delta t}{\tau}\right)f_{i}^{\rm Reg}+\Delta t{F}_{i}.$ (17)

###  A.2 Kinetic Boundary Condition

The no-slip boundary conditions based on bounce back mechanism, where the
directions of the incoming distribution functions are simply reversed when it
encounters the boundary node, as the name suggests fails as to uncover the
slip velocity at the boundary at finite Kn. To overcome this shortcoming of
the method, a diffusively reflecting solid wall boundary condition was
introduced by Ansumali and Karlin ansumali2002kinetic . Based on kinetic
theory interpretation, this boundary condition redistributes the population
coming towards the wall in such a way that mass-balance and normal-flux
conditions are fulfilled. In the discrete sense, the distribution function at
the boundary/wall takes the following form

$f_{i}({\boldsymbol{x}}_{w},t)=Kf_{i}^{\rm eq}(\rho,{\boldsymbol{u}}_{w}),$
(18)

where subscript $w$ denotes the wall and

$K=\frac{\sum_{{\boldsymbol{c}}_{i}\cdot{\boldsymbol{n}}<0}|({\boldsymbol{c}}_{i}-{\boldsymbol{u}}_{w})\cdot{\boldsymbol{n}}|f_{i}}{\sum_{{\boldsymbol{c}}_{i}\cdot{\boldsymbol{n}}>0}|({\boldsymbol{c}}_{i}-{\boldsymbol{u}}_{w})\cdot{\boldsymbol{n}}|f_{i}^{\rm
eq}(\rho,{\boldsymbol{u}}_{w})},$ (19)

with ${\boldsymbol{n}}$ being the unit normal direction.

(a) $n_{x}$

(b) $n_{y}$

Figure 9: Unit normals at the solid boundaries

The term K can be understood as the ratio of outgoing flux from the wall and
incoming equilibrium flux coming towards the wall. When the boundaries are
stationary (${\boldsymbol{u}}_{w}=0$), the term K can be written as:

$K=\frac{\sum_{A_{i}<0}|A_{i}|f_{i}}{\sum_{A_{i}>0}|A_{i}|f_{i}^{\rm eq}},$
(20)

with $A_{i}={c}_{ix}{n}_{x}+{c}_{iy}{n}_{y}$. To implement this boundary
condition one needs the information of unit normals at the boundaries. Fig
9(a) and 9(a) respectively represents the unit normals in the $x$-direction,
$n_{x}$, and the same in $y$-direction, $n_{y}$ singh2017impact ;
singh2017influence  for a circular obstacle. In the present study, the normals
for an individual circle with center $(x_{\rm cen},y_{\rm cen})$, are
calculated as

$\hat{\boldsymbol{n}}=\frac{{\boldsymbol{x}}-{\boldsymbol{x}}_{\rm
cen}}{|{\boldsymbol{x}}-{\boldsymbol{x}}_{\rm cen}|}.$ (21)

## References

  * [1] J Bear.  Dynamics of fluids in porous media.  1972\.
  * [2] Behzad Ghanbarian, Allen G Hunt, Robert P Ewing, and Muhammad Sahimi.  Tortuosity in porous media: a critical review.  Soil science society of America journal, 77(5):1461–1477, 2013\.
  * [3] Bert Metz, Ogunlade Davidson, HC De Coninck, Manuela Loos, and Leo Meyer.  IPCC special report on carbon dioxide capture and storage.  Cambridge: Cambridge University Press, 2005.
  * [4] Larry W Lake.  Enhanced oil recovery.  1989\.
  * [5] Philip Crosbie Carman.  Fluid flow through a granular bed.  Trans. Inst. Chem. Eng. London, 15:150–156, 1937.
  * [6] Francis AL Dullien.  Porous media: fluid transport and pore structure.  Academic press, 2012.
  * [7] Antti Koponen, M Kataja, and JV Timonen.  Tortuous flow in porous media.  Physical Review E, 54(1):406, 1996.
  * [8] Maciej Matyka, Arzhang Khalili, and Zbigniew Koza.  Tortuosity-porosity relation in porous media flow.  Physical Review E, 78(2):026306, 2008.
  * [9] Artur Duda, Zbigniew Koza, and Maciej Matyka.  Hydraulic tortuosity in arbitrary porous media flow.  Physical Review E, 84(3):036319, 2011.
  * [10] Farzam Javadpour, Douglas Fisher, and Michelle Unsworth.  Nanoscale gas flow in shale gas sediments.  Journal of Canadian Petroleum Technology, 46(10), 2007.
  * [11] CM Freeman, GJ Moridis, and TA Blasingame.  A numerical study of microscale flow behavior in tight gas and shale gas reservoir systems.  Transport in porous media, 90:253–268, 2011.
  * [12] LJ Klinkenberg.  The permeability of porous media to liquids and gases: American petroleum institute drilling and production practice paper 41-200: Washington.  DC, American Petroleum Institute, 1941.
  * [13] Roberto Benzi, Sauro Succi, and Massimo Vergassola.  The lattice boltzmann equation: theory and applications.  Physics Reports, 222(3):145–197, 1992.
  * [14] Shiyi Chen and Gary D Doolen.  Lattice boltzmann method for fluid flows.  Annual review of fluid mechanics, 30(1):329–364, 1998.
  * [15] Sauro Succi.  The lattice Boltzmann equation: for fluid dynamics and beyond.  Oxford university press, 2001.
  * [16] Cyrus K Aidun and Jonathan R Clausen.  Lattice-boltzmann method for complex flows.  Annual review of fluid mechanics, 42:439–472, 2010.
  * [17] T Krüger et al.  The lattice boltzmann method: Principles and practice. sl: Springer international publishing.  2017\.
  * [18] Y. Zhou, R. Zhang, I. Staroselsky, H. Chen, W.T. Kim, and M.S. Jhon.  Simulation of micro-and nano-scale flows via the lattice Boltzmann method.  Physica A: Statistical Mechanics and its Applications, 362(1):68–77, 2006.
  * [19] R. Zhang, X. Shan, and H. Chen.  Efficient kinetic method for fluid simulation beyond the Navier-Stokes equation.  Physical Review E, 74(4):046703, 2006.
  * [20] J. Latt and B. Chopard.  Lattice Boltzmann method with regularized pre-collision distribution functions.  Mathematics and Computers in Simulation, 72(2):165–168, 2006.
  * [21] Shi Tao and Zhaoli Guo.  Boundary condition for lattice boltzmann modeling of microscale gas flows with curved walls in the slip regime.  Physical Review E, 91(4):043305, 2015.
  * [22] Zhaoli Guo, Jishun Qin, and Chuguang Zheng.  Generalized second-order slip boundary condition for nonequilibrium gas flows.  Physical Review E, 89(1):013021, 2014.
  * [23] Santosh Ansumali and Iliya V Karlin.  Kinetic boundary conditions in the lattice boltzmann method.  Physical Review E, 66(2):026311, 2002.
  * [24] A Montessori, P Prestininzi, M La Rocca, and S Succi.  Lattice boltzmann approach for complex nonequilibrium flows.  Physical Review E, 92(4):043308, 2015.
  * [25] Shiwani Singh, Fei Jiang, and Takeshi Tsuji.  Impact of the kinetic boundary condition on porous media flow in the lattice boltzmann formulation.  Physical Review E, 96(1):013303, 2017.
  * [26] Shiwani Singh, Fei Jiang, and Takeshi Tsuji.  Influence of slip flow at fluid-solid interface upon permeability of natural rock.  Energy Procedia, 114:3572–3577, 2017.
  * [27] MA Knackstedt and X Zhang.  Direct evaluation of length scales and structural parameters associated with flow in porous media.  Physical Review E, 50(3):2134, 1994.
  * [28] Xiaodong Zhang and Mark A Knackstedt.  Direct simulation of electrical and hydraulic tortuosity in porous solids.  Geophysical research letters, 22(17):2333–2336, 1995.
  * [29] Matteo Icardi, Gianluca Boccardo, Daniele L Marchisio, Tiziana Tosco, and Rajandrea Sethi.  Pore-scale simulation of fluid flow and solute dispersion in three-dimensional porous media.  Physical review E, 90(1):013032, 2014.
  * [30] Maciej Matyka, Jarosław Gołembiewski, and Zbigniew Koza.  Power-exponential velocity distributions in disordered porous media.  Physical Review E, 93(1):013110, 2016.
  * [31] Bagus P Muljadi, Martin J Blunt, Ali Q Raeini, and Branko Bijeljic.  The impact of porous media heterogeneity on non-darcy flow behaviour from pore-scale simulation.  Advances in water resources, 95:329–340, 2016.
  * [32] AA El-Zehairy, M Mousavi Nezhad, Vahid Joekar-Niasar, I Guymer, N Kourra, and MA Williams.  Pore-network modelling of non-darcy flow through heterogeneous porous media.  Advances in Water Resources, 131:103378, 2019.
  * [33] Krzysztof M Graczyk and Maciej Matyka.  Predicting porosity, permeability, and tortuosity of porous media from images by deep learning.  Scientific reports, 10(1):21488, 2020.
  * [34] Maciej Matyka and Zbigniew Koza.  How to calculate tortuosity easily?  In AIP Conference Proceedings 4, volume 1453, pages 17–22. American Institute of Physics, 2012.
  * [35] S Mohammed and T Reis.  Lattice boltzmann method with moment-based boundary conditions for rarefied flow in the slip regime.  Physical Review E, 104(4):045309, 2021.
  * [36] Ali Beskok and George Em Karniadakis.  Report: a model for flows in channels, pipes, and ducts at micro and nano scales.  Microscale thermophysical engineering, 3(1):43–77, 1999.
  * [37] Faruk Civan.  Effective correlation of apparent gas permeability in tight porous media.  Transport in porous media, 82:375–384, 2010.
  * [38] Lei Wu, Haihu Liu, Jason M Reese, and Yonghao Zhang.  Non-equilibrium dynamics of dense gas under tight confinement.  Journal of Fluid Mechanics, 794:252–266, 2016.
  * [39] Lei Wu, Minh Tuan Ho, Lefki Germanou, Xiao-Jun Gu, Chang Liu, Kun Xu, and Yonghao Zhang.  On the apparent permeability of porous media in rarefied gas flows.  Journal of Fluid Mechanics, 822:398–417, 2017.
  * [40] Guang Yang and Bernhard Weigand.  Investigation of the klinkenberg effect in a micro/nanoporous medium by direct simulation monte carlo method.  Physical Review Fluids, 3(4):044201, 2018.
  * [41] Wei Su, Minh Tuan Ho, Yonghao Zhang, and Lei Wu.  Gsis: An efficient and accurate numerical method to obtain the apparent gas permeability of porous media.  Computers & Fluids, 206:104576, 2020.
  * [42] Felix J Meigel, Thomas Darwent, Leonie Bastin, Lucas Goehring, and Karen Alim.  Dispersive transport dynamics in porous media emerge from local correlations.  Nature Communications, 13(1):5885, 2022.
  * [43] Fei Jiang and Takeshi Tsuji.  Changes in pore geometry and relative permeability caused by carbonate precipitation in porous media.  Physical Review E, 90(5):053306, 2014.
  * [44] Fei Jiang, Yaotian Guo, Takeshi Tsuji, Yoshitake Kato, Mai Shimokawara, Lionel Esteban, Mojtaba Seyyedi, Marina Pervukhina, Maxim Lebedev, and Ryuta Kitamura.  Upscaling permeability using multiscale x-ray-ct images with digital rock modeling and deep learning techniques.  Water Resources Research, 59(3):e2022WR033267, 2023.
  * [45] Prabhu Lal Bhatnagar, Eugene P Gross, and Max Krook.  A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems.  Physical review, 94(3):511, 1954.
  * [46] Zhaoli Guo, Chuguang Zheng, and Baochang Shi.  Discrete lattice effects on the forcing term in the lattice boltzmann method.  Physical review E, 65(4):046308, 2002.


# Extraction of kinetics from equilibrium distributions of states using the
Metropolis Monte Carlo method

Sergei F. Chekmarev  Institute of Thermophysics, SB RAS, 630090 Novosibirsk,
Russia  Department of Physics, Novosibirsk State University, 630090
Novosibirsk, Russia<EMAIL_ADDRESS>
###### Abstract

The Metropolis Monte Carlo (MC) method is used to extract reaction kinetics
from a given equilibrium distribution of states of a complex system. The
approach is illustrated by the folding/unfolding reaction for two proteins - a
model $\beta$-hairpin and $\alpha$-helical protein $\alpha_{3}$D. For
$\beta$-hairpin, the free energy surfaces (FESs) and free energy profiles
(FEPs) are employed as the equilibrium distributions of states, playing a role
of the potentials of mean force to determine the acceptance probabilities of
new states in the MC simulations. Based on the FESs and PESs for a set of
temperatures that were simulated with the molecular dynamics (MD) method, the
MC simulations are performed to extract folding/unfolding rates. It has been
found that the rate constants and first-passage time (FPT) distributions
obtained in the MC simulations change with temperature in good agreement with
those from the MD simulations. For $\alpha_{3}$D protein, whose equilibrium
folding/unfolding was studied with the single-molecule FRET method (Chung et
al., J. Phys. Chem. A, 115, 2011, 3642), the FRET-efficiency histograms at
different denaturant concentrations were used as the equilibrium distribution
of protein states. It has been found that the rate constants for folding and
unfolding obtained in the MC simulations change with denaturant concentration
in reasonable agreement with the constants that were extracted from the photon
trajectories on the basis of theoretical models. The simulated FPT
distributions are single-exponential, which is consistent with the assumption
of two-state kinetics that was made in the theoretical models. The promising
feature of the present approach is that it is based solely on the equilibrium
distributions of states, without introducing any additional parameters to
perform simulations, which suggests its applicability to other complex
systems.

††preprint: AIP/123-QED

##  I Introduction

The equilibrium distribution of states and the kinetics of their
interconversion are two inherent characteristics of any complex system that
are important to know for understanding its behavior. At the same time, the
former are generally much easier to determine than the latter, both in
experiment and in simulations. One representative example is the protein
folding/unfolding reaction. Single-molecule experiments, such as the Förster
resonance energy transfer (FRET), atomic-force microscopy (AFM) and optical
tweezers (OT) techniques, make it possible to determine the equilibrium
distributions of states, either in the form of the FRET-efficiency histograms
(FRET schuler2002probing ; huang2009direct ; chung2011extracting ;
liu2012exploring ; schuler2013single ; eaton2021modern , folding/unfolding) or
free energy profiles (AFM yu2012energy ; batchelor2020protein  and OTT
bustamante2020single , unfolding). However, in order to extract kinetics from
these distributions, theoretical models are needed (FRET mckinney2006analysis
; gopich2009decoding ; gopich2010fret ; ramanathan2015method ;
osadko2016amplitude , and AFM and OTT hummer2001free ; hummer2005free ;
hyeon2007measuring ; dudko2008theory ). A similar difference in the
possibility of obtaining the equilibrium distribution of states and
determining the kinetics exists for the molecular dynamics (MD) simulations of
folding of relatively large proteins (of the order of one hundred of residues
and larger). While brute force all-atom simulations, which are able to give a
detailed insight into the folding kinetics, are very time-consuming
shaw2010atomic ; lindorff-larsen2011how , several methods exist that allow
very efficient calculation of the equilibrium distributions of protein states
camilloni2018advanced , such as the replica exchange molecular dynamics
sugita1999replica ; rhee2003multiplexed ; rao2003replica  and the weighted
histogram analysis method ferrenberg1989  in combination with sampling from
biased molecular dynamics simulations shea2001from . Therefore, it is tempting
to find out whether the knowledge of the equilibrium distribution of states of
a system allows us to obtain information about its kinetics, and if so, to
what extent.

The present paper examines the possibility of extracting kinetics from the
equilibrium distributions of protein states with the use of the Metropolis’s
Monte Carlo (MC) algorithm metropolis1953equation ; frenkel2004introduction ,
assuming that the acceptance probability of a new state can be determined by
the equilibrium state distribution known from simulations or experiment. The
typical examples of such distributions are free energy surfaces (FESs) and
free energy profiles (FEPs), although, in general, any equilibrium
distribution of states along one or several collective variables can be used
for this purpose. In contrast to the nonequilibrium distributions of protein
states that determine the probability fluxes of folding
chekmarev2008hydrodynamic ; kalgin2014first , the equilibrium distributions
represent potentials of mean-force (PMF) roux1995calculation ;
trzesniak2007comparison ; shea2001from , and thus they can be assumed to
govern kinetics. A known shortcoming of the MC methods is that a generated MC
trajectory represents a Markovian sequence of states rather than the change of
these states in time. Therefore, only the relative values of the reaction
rates for different reaction conditions, e.g., for temperature, are available.
In order to have absolute value of rates, a fit to the time-scale of the
original data is required huitema1999can ; rutkai2010dynamic ; sanz2010dynamic
; bal2014time . At the same time, in contrast to other possible approaches
which could employ the PMFs to calculate kinetics, e.g., the use of Fokker-
Planck or Langevin equations kikuchi1991metropolis ; tiana2007use , no
additional parameters are required to perform the MC simulations for a given
equilibrium distribution of states. It is a considerable advantage of the
present approach.

The approach is illustrated for two systems. One is a model $\beta$-hairpin
protein. In order to have baseline results, the equilibrium MD simulations
were performed. The obtained FESs and FEPs were employed as the input
distributions of states for subsequent MC simulations. It has been shown that
the rate constants of folding and unfolding calculated from the FESs and FEPs
change with temperature in good agreement with those from the MD simulations.
The other system is an $\alpha$-helical protein $\alpha_{3}$D, for which
folding and unfolding rate constants for different denaturant concentrations
were extracted from single-molecule FRET trajectories chung2011extracting
using theoretical models gopich2009decoding ; gopich2010fret . In this case,
the FRET-efficiency histograms (FEHs) obtained in the experiment were employed
as the equilibrium distributions of states for the MC simulations. It has been
found that the folding and unfolding rate constants obtained in the MC
simulations vary with the denaturant concentration in reasonable agreement
with those extracted from the FRET trajectories. For both proteins, the first-
passage time (FPT) distributions in the MC simulations have been found to be
single-exponential, which agree with those obtained in the MD simulations
($\beta$-hairpin protein) and correspond to the assumption of two-state
kinetics employed in the theoretical models to extract rate constants
($\alpha_{3}$D protein).

The paper is organized as follows. Section II describes the MD (Sect. II.1)
and MC (Sect. II.2) simulations for the $\beta$-hairpin protein, Sect. III
present the results of the MC simulations for the $\alpha_{3}$D protein, and
Sect. IV briefly summarizes the results of the work.

##  II  $\beta$-hairpin Protein: System and Simulation

###  II.1 System and Molecular Dynamics Simulations

To perform MD simulations, a coarse-grained (Cα-bead) protein model and a Gō-
type interaction potential go1983theoretical  were used. Briefly, the approach
is as follows (for details, see Ref. chekmarev2013protein, ). The Cα-model was
constructed on the basis of the solution NMR data for 12-residue HP7 protein
(2evq.pdb) andersen2006minimization . Specifically, the first conformation in
the NMR ensemble of the protein structures was employed for this purpose,
which was then considered as a reference native structure of the protein. The
Gō-type potential accounted for the rigidity of the backbone and the
contributions of native and non-native contacts hoang2000molecular . Two Cα-
beads were considered to be in native contact if they were not the nearest
neighbors along the protein chain and had the interbead distance not longer
than $d_{\mathrm{cut}}=7.5$Å. In this case, the number of native contacts is
$N_{\mathrm{nat}}=N_{\mathrm{nat}}^{\mathrm{NAT}}=27$. The simulations were
performed with a constant-temperature MD based on the coupled set of Langevin
equations biswas1986simulated . The time-step was $\Delta t=0.0125\tau$, where
$\tau$ is the characteristic time. At the length scale $l=7.5\mathrm{\AA}$Å
and the attractive energy $\epsilon=2.2$ kcal/mol miyazawa1996residue ,
$\tau=(Ml^{2}/\epsilon)^{1/2}\approx 2.6$ ps, where $M=110$ Da is the average
mass of the residue. The friction constant $\gamma$ in the Langevin equations
was $\gamma=10M/\tau$. In what follows, the length, in particular, the radius
of gyration, is measured in angstroms. Unless otherwise noted, the other
quantities (energy, temperature and time) are dimensionless; specifically, the
energy and temperature are measured in units $\epsilon$ (in the latter case
with the Boltzmann constant set to unity), and the time is measured in units
of $\tau$.

The equilibrium simulations were carried out for a range of temperatures at
which both folding and unfolding events occurred quite frequently,
specifically, at $T=$ 0.35, 0,375, 0.4, 0.425, and 0.45. The regime of folding
varied from downhill (low temperatures) to uphill (larger temperatures)
folding through a marginal two-state folding (intermediate temperatures). For
each temperature, the MD trajectory was run until $10^{3}$ folding/unfolding
events took place. The protein was considered to be unfolded if
$N_{\mathrm{nat}}\leq 5$, and to be folded if the root-mean-square-deviation
(RMSD) from the reference native structure was less than 1Å chekmarev2021first
. Based on the simulated trajectories, the FESs and FEPs were calculated. The
FES was constructed using the number of native contacts ($N_{\mathrm{nat}}$)
and the radius of gyration ($R_{\mathrm{g}}$) as collective variables
best2013native ; shea2001from

$F(N_{\mathrm{nat}},R_{\mathrm{g}})=-T\ln P(N_{\mathrm{nat}},R_{\mathrm{g}})$
(1)

where $P(N_{\mathrm{nat}},R_{\mathrm{g}})$ is the probability to find the
protein in the point $(N_{\mathrm{nat}},R_{\mathrm{g}})$. The FEP was
calculated as

$F(N_{\mathrm{nat}})=-T\ln P(N_{\mathrm{nat}})$ (2)

where $P(N_{\mathrm{nat}})$ is the probability that the protein has
$N_{\mathrm{nat}}$ native contacts. The resulting FESs and FEPs for the lower,
intermediate and upper temperatures are shown in Fig. 1a-c and Fig. 2,
respectively. In structure, they are common for folding of $\beta$-hairpins
(experiment munoz1997folding , and simulation studies in explicit zhou2001free
; gao2015correct  and implicit dinner1999understanding ; zagrovic2001beta
solvents).

Figure 1: $\beta$-hairpin, free energy surfaces: (a) $T=0.35$, (b) $T=0.4$,
and (c) $T=0.45$. The radius of gyration is given in angstroms.

Figure 2: $\beta$-hairpin, free energy profiles: $T=0.35$ (blue and cyan),
$T=0.4$ (black and gray), and $T=0.45$ (red and magenta). The first colors in
the brackets are for the MD simulations, and second colors are for the MC
simulations based on the FESs of Fig. 1.

Both folding and unfolding reactions were considered. For this, the
equilibrium MD trajectory was divided into segments of the trajectory
corresponding to the transition from the unfolded to the native state
(folding) and to the backward transition (unfolding) chekmarev2021first .
Based on these segments of the trajectory, the MD rate constants of folding
and unfolding, as well as the corresponding FPT distributions, were calculated
(see below).

###  II.2  Monte Carlo Simulations

Having the FESs and FEPs, the MC simulations were performed. According to the
Metropolis algorithm metropolis1953equation ; frenkel2004introduction , the
acceptance probability of the transition from state A to B is determined as

$W(\mathrm{A}\rightarrow\mathrm{B})=\begin{cases}1&{\text{if}}\,\,\,\,P(\mathrm{B}/P(\mathrm{A})\geq
1\\\
P(\mathrm{B})/P(\mathrm{A})&{\text{if}}\,\,\,\,P(\mathrm{B})/P(\mathrm{A})<1\end{cases}$
(3)

where $P(\mathrm{A})$ and $P(\mathrm{B})$ are the equilibrium probabilities of
the old ($\mathrm{A}$) and new ($\mathrm{B}$) states, respectively. A trial
move from $\mathrm{A}$ to $\mathrm{B}$ is accepted if
$P(\mathrm{B})/P(\mathrm{A})\geq\alpha$, where $\alpha$ is a random number
uniformly distributed between 0 and 1. In application to the free energy
landscapes of Figs. 1 and 2,
$P(\mathrm{B})/P(\mathrm{A})=\exp\\{-[F(\mathrm{B})-F(\mathrm{A})]/T\\}$. For
FESs (Fig. 1), $\mathrm{A}$ and $\mathrm{B}$ represent the points in the
$(N_{\mathrm{nat}},R_{\mathrm{g}})$ space, and for FEPs (Fig. 2), the points
in the $N_{\mathrm{nat}}$ space. The corresponding free energies
$F(\mathrm{A})$ and $F(\mathrm{B})$ are determined by Eqs. (1) and (2).

In the case of FESs, the MC simulation space represented a rectangle with
$0\leq N_{\mathrm{nat}}\leq 28$ and $0\leq R_{\mathrm{g}}\leq 15$Å. It was
divided into $28\times 30$ segments ($dN_{\mathrm{nat}}=1$ and
$dR_{\mathrm{g}}=0.5$), forming a mesh of $29\times 31$ discrete states. The
trial moves consisted of random shifts of the point along the
$N_{\mathrm{nat}}$ and $R_{\mathrm{g}}$ axes by -1 or +1 mesh step so that the
8 nearest points surrounding the current point were sampled with equal
probability. To verify that detailed balance is satisfied, a long
("equilibrium") MC trajectory with many folding/unfolding events was run. It
was found that the balance was approximately maintained for each point-to-
point transition and was exactly fulfilled for the transitions within the
simulation cell, i.e. the number of transitions from the central point of the
cell to the eight surrounding points was equal to the number of the backward
transitions. Accordingly, the equilibrium MC simulations successfully
reproduced the original equilibrium distributions on which they were based
(one example is shown in Fig. 2). To simulate the process of folding, the
trajectories were initiated at point (5, 11.8), representing the unfolded
state, and terminated at point (27, 5.78), representing the native state,
(Fig. 1a-c). For the unfolding trajectories, these points served as the
terminal and initial states, respectively.

In the case of FEPs, the MC simulations were conducted in the linear
$N_{\mathrm{nat}}$ space with the boundaries $N_{\mathrm{nat}}=0$ and
$N_{\mathrm{nat}}=27$. The trial moves represented random changes of the
current value of $N_{\mathrm{nat}}$ by -1 or +1, which were taken with equal
probability. In equilibrium, detailed balance was strictly maintained. The
folding trajectories started at $N_{\mathrm{nat}}=5$ and terminated at
$N_{\mathrm{nat}}=27$, and the unfolding trajectories used these points as
terminal and initial state, respectively.

Figure 3: $\beta$-hairpin, rate constants of (a) folding and (b) unfolding:
the MD simulations (black diamonds and dashed line), the FES based MC
simulations (read triangles and dotted line), and the FEP based MC simulations
(blue triangles and dashed line). The lines are to guide the eye.

Figure 4: $\beta$-hairpin, first-passage time distributions: (a) folding, and
(b) unfolding. Solid lines are for the MD simulations, dashes lines are for
the MC simulations based on the free energy surfaces of Fig. 1, and dotted
lines are for the MC simulations based on the free energy profiles of Fig. 2.
Temperatures and lines: $T=0.35$ (blue), $T=0.375$ (olive), $T=0.3$ (black),
$T=0.425$ (magenta), and $T=0.45$ (red).

In each case (the FES or FEP based simulations, the protein folding or
unfolding, and a given temperature), $10^{3}$ trajectories were run. For the
FES based simulations, the acceptance probability increased from $\approx
0.57$ ($T=0.35$) to $\approx 0.67$ ($T=0.45$) for folding trajectories, and
from $\approx 0.57$ ($T=0.35$) to $\approx 0.62$ ($T=0.45$) for unfolding
trajectories, and for the FEP based simulations, it increased from $\approx
0.87$ ($T=0.35$) to $\approx 0.9$ ($T=0.45$) for both folding and unfolding
trajectories.

In order to fit the MC steps to the MD time steps, the average ratios btween
the MC and MD rate constants over the whole temperature range were calculated.
In particular, for the FES based folding trajectories, the time step was
determined as $(\Delta
t)^{\mathrm{F}}_{\mathrm{FES}}=(1/N)\sum_{i=1}^{i=N}(k_{\mathrm{F}})_{\mathrm{MC},i}/(k_{\mathrm{F}})_{\mathrm{MD},i}$,
where $N=5$ is the number of temperatures, and
$(k_{\mathrm{F}})_{\mathrm{MC}}$ and $(k_{\mathrm{F}})_{\mathrm{MD}}$ are the
rate constants for the MC and MD folding trajectories, respectively. The rate
constants were calculated as the inverses of the corresponding mean first-
passage times (MFPTs) in units of number of steps (the MC simulations) and
$\tau$ (the MD simulations). Similarly, the time steps for the FEP based and
unfolding MC simulations were determined. It was found that $(\Delta
t)^{\mathrm{F}}_{\mathrm{FES}}\approx 3.4\times 10^{-3}\tau$, $(\Delta
t)^{\mathrm{F}}_{\mathrm{FEP}}\approx 6.3\times 10^{-3}\tau$, and $(\Delta
t)^{\mathrm{U}}_{\mathrm{FES}}\approx(\Delta
t)^{\mathrm{U}}_{\mathrm{FEP}}\approx 0.18\tau$. These figures show that the
MC time scales generally depend on both the direction of the process (folding
or unfolding) and the dimensionality of the state distributions for the MC
simulations (FES or FEP).

Figures 3 and 4 compare the results of the MC and MD simulations. Figure 3
presents the rate constants obtained in the MC simulations in comparison to
those in the MD simulations. The MC "time” scales were adjusted to the MD time
scales using the above values of the MC steps in $\tau$ units, e.g., for the
MC folding simulations based on the FES, the adjusted rate constant at
$T=0.35$ is obtained to be
$(k_{\mathrm{F}})_{\mathrm{MC}}=(\tilde{k}_{\mathrm{F}})_{\mathrm{MC}}/(\Delta
t)^{\mathrm{F}}_{\mathrm{FES}}\approx 2.6\times 10^{-3}$, where
$(\tilde{k}_{\mathrm{F}})_{\mathrm{MC}}\approx 8.75\times 10^{-6}$ is the
original rate constant in the inverse number of the MC steps, and $(\Delta
t)^{\mathrm{F}}_{\mathrm{FES}}\approx 3.4\times 10^{-3}\tau$. Figure 3 shows
that the rate constants from the MC simulations are in good agreement with
those from the MD simulations through all temperature range. The same result
would evidently be obtained if the MC step in $\tau$ units was determined by
equating the MC rate constant at one temperature with the corresponding MD
rate constant. The above values of the time steps were also used to adjust the
time scales for the MC FPT distributions. Figure 4 presents these
distributions in comparison with the distributions obtained in the MD
simulations and shows that, similar to the rate constants, they are also in
agreement, including a single-exponential decay. Importantly, the good
agreement of the MC rate constants with the MD constants through all
temperature range evidences that the time step in the MC simulations is not
affected by temperature.

##  III $\alpha_{3}$D Protein: Experiment and Monte Carlo Simulations

The equilibrium folding/unfolding of this 73-residue three-helix bundle
protein (2a3d.pdb) at room temperature and different denaturant (GdmCl)
concentrations (1.5 M to 3 M) was studied using single-molecule FRET
chung2011extracting . Two models assuming two-state kinetics were employed to
analyze the photon trajectories. One is a Maximum Likelihood Estimation (MLE),
in which the photon trajectory were fitted to the sequences of transitions
between unfolded and folded states gopich2009decoding , and the other is an
approximation of the FRET-efficiency histogram (FEH) by a sum of Gaussians
whose weights depended on the probabilities of surviving the folded and
unfolded states for the bin time (3G) gopich2010fret . Using these approaches,
the rate constants of folding and unfolding were determined, which were found
in good agreement, except for low values of the denaturant concentration,
specifically for 1.5 M and, partly, for 2.0 M (cf. the corresponding lines of
Table 1 of Ref. chung2011extracting, ). Figure 5 shows the FEHs for freely

Figure 5: $\alpha_{3}$D: FRET-efficiency histograms chung2011extracting .
Figures inside the panels represent the GdmCl concentrations. Labels U and F
indicate the unfoded and folded state, respectively.

diffusing protein molecules. The time bin to count donor and acceptor photons
is 2 ms. Two peaks are significant - at intermediate and higher values of
FRET-efficiency ($E$), which represent the unfolded (U) and folded (F) states,
as indicated in Fig. 5. The peaks at $E\sim 0$ originate from molecules with
inactive acceptors and are not included in consideration. It was found that
the FEHs are well approximated by a sum of three Gaussians (3G model). Two
Gaussians describe the unfolded and folded states, and the third one accounts
for the intermediate values of $E$ that appear due to the transitions between
the folded and unfolded states.

Figure 6: $\alpha_{3}$D: Comparison of the rate constants for protein (a)
folding and (b) unfolding. Open triangles show the experimental rate constants
obtained with the 3G model chung2011extracting , and solid triangles represent
simulated rate constants.

Figure 7: $\alpha_{3}$D: Simulated first-passage time distributions for
protein (a) folding and (b) unfolding.

It should be noted that the present FEHs may not represent the true PMFs, as
it would be if they obeyed the Förster "distance-efficiency" equation
$E=1/[1+(r/R_{0})^{6}]$ ($r$ is the donor-to-acceptor separation distance, and
$R_{0}$ is the Förster radius). First, the number of photons in time bin
($\sim 10^{2}$) is not sufficiently large to avoid the effect of shot noise on
the distribution of protein states. In this case, the broadening of the folded
and unfolded peaks may mainly be a results of short noise rather than of the
presence of folded- and unfolded-like conformations of the protein. Second,
the binning time (2 ms) is comparable to the relaxation time ($\sim 1$ms), so
that the signal from different protein conformations, which typically
interconvert on a much shorter time scale ($\sim 1\mu$s), may be placed in the
same bin. In particular, the intermediate values of $E$ between the unfolded
and folded peaks may represent a mixture of photons that come from the
unfolded and folded states rather than the transitional conformations.
However, if detailed balance is fulfilled, it is more important that the
intermediate values of $E$ adequately represent the probabilities of the
intermediate protein conformations rather than the conformations themselves.

The MC moves were determined by Eq. (3), in which the values of FRET-
efficiency represented the protein states, and the numbers of bursts in the
FEHs served as (non-normalized) probabilities of these states, $P(E)$.
According to the FEHs of Fig. 5, the whole range of FRET-efficiencies, from 0
to 1, was divided into 31 segments. To take into account that the peak at
$E\sim 0$ corresponds to the inactive acceptor, the populations of states at
$0<E<E_{\mathrm{thr}}$, where $E_{\mathrm{thr}}$ is the value of $E$ at which
$P(E)$ has a minimum between the inactive and unfolded peaks, were set to
zero, $P(E)=0$. Each MC step was taken to be an equally probable random move
to a neighboring segment. Detailed balance was strictly satisfied similar to
the MC simulations for $\beta$-hairpin (Sect. II.2). The unfolded and folded
states were chosen to correspond to the maxima of the unfolded and folded
peaks (Fig. 5), with the values of FRET-efficiency $E_{\mathrm{U}}=0.61$ and
$E_{\mathrm{F}}=0.93$, respectively. To simulate folding, the MC trajectories
were initiated at $E=E_{\mathrm{U}}$ and terminated at $E=E_{\mathrm{F}}$, and
to simulate unfolding, the states $E_{\mathrm{F}}$ and $E_{\mathrm{U}}$ were
used as the initial and terminal states, respectively. In all cases, i.e., at
each denaturant concentration for folding and unfolding, $10^{3}$ MC
trajectories were run. For both folding and unfolding reactions, the
acceptance probability varied insignificantly, from $\approx 0.87$ to $\approx
0.92$. Based on the simulated folding and unfolding trajectories, the rate
constants for folding ($k_{\mathrm{F}}$) and unfolding ($k_{\mathrm{U}}$) were
determined, and the corresponding FPT distributions were constructed. Similar
to $\beta$-hairpin, the MC “time” scale was adjusted to the experimental time
scale by calculating the average ratios between simulated and experimental
rate constants. Specifically, the MC time steps for folding and unfolding were
determined, respectively, as $(\Delta
t)^{\mathrm{F}}=(1/N)\sum_{i=1}^{i=N}(k_{\mathrm{F}})_{\mathrm{MC},i}/(k_{\mathrm{F}})_{\mathrm{exp},i}$
and $(\Delta
t)^{\mathrm{U}}=(1/N)\sum_{i=1}^{i=N}(k_{\mathrm{U}})_{\mathrm{MC},i}/(k_{\mathrm{U}})_{\mathrm{exp},i}$,
where $N=3$ is the number of denaturant concentrations ranging from 2.0 M to
2.5 M (i.e., 2.0 M, 2.25 M, and 2.5 M), at which the relative populations of
both folded and unfolded states were larger than $1\%$. It was obtained
$(\Delta t)^{\mathrm{F}}\approx 7.5\times 10^{-3}$ms and $(\Delta
t)^{\mathrm{U}}\approx 8.7\times 10^{-3}$ms. Given the MC time steps, the rate
constants for folding and unfolding were recalculated as
$(k_{\mathrm{F}})_{\mathrm{MC}}=(\tilde{k}_{\mathrm{F}})_{\mathrm{MC}}/(\Delta
t)^{\mathrm{F}}$ and
$(k_{\mathrm{U}})_{\mathrm{MC}}=(\tilde{k}_{\mathrm{U}})_{\mathrm{MC}}/(\Delta
t)^{\mathrm{U}}$, where the rate constants denoted by tilde are the original
MC constants in the inverse number of the MC steps.

Figure 6a-b compares the adjusted MC rate constants with those extracted from
the FRET-trajectories with the 3G model (Table 1 of Ref. chung2011extracting,
). It shows that the MC simulations based on the FRET-histograms reasonably
predict the change of the rate constants with denaturant concentration, for
both folding and unfolding reactions. The only exception are the
concentrations at which the normalized population of the target state is as
small as $\approx 1\%$, as at 3M GdmCl in the case of folding, where
$P(0.93)\approx 1.2\times 10^{-2}$, and at 1.5M GdmCl in the case of
unfolding, where $P(0.61)\approx 9\times 10^{-3}$. More generally, the MC
simulations are expected to give satisfactory results when the peaks for the
target states are resolved. Figure 7a-b also presents the simulated FPT
distributions for folding and unfolding. The distributions are single-
exponential, which is consistent with the assumption of two-state kinetics
that was made in the theoretical models for extracting the rate constants from
the FRET data chung2011extracting .

##  IV Conclusions

It has been shown that the Metropolis MC simulations based on the equilibrium
distributions of states can be successfully used to predict how characteristic
reaction times in a complex system change with the determining conditions. As
an example, the process of folding and unfolding for two proteins was
considered. The first protein, a $\beta$-hairpin protein, whose coarse-grained
(Cα-bead) model was constructed on the basis of the solution NMR data for
12-residue HP7 protein andersen2006minimization , was employed to develop and
test the approach. Performing MD simulations, the baseline characteristics of
the folding/unfolding process were determined for a set of temperatures - the
equilibrium free energy surfaces (FESs) and free energy profiles (FEPs), the
rate constants of folding and unfolding, and the first-passage time (FPT)
distributions for these reactions. Using the FESs and FEPs as the potentials
of mean-force (PMFs), the MC simulations were carried out to obtain relative
rate constants of folding and unfolding and the FPT distributions. It has been
found that they change with temperature in good agreement with those from the
MD simulations and, after fitting the MC steps to the MD time steps, they
agree well in absolute value. The other protein is a three-helical
$\alpha_{3}$D protein, whose folding/unfolding reaction was studied at
different different denaturant concentrations using single-molecule FRET
method chung2011extracting . In this case, the measured FRET-efficiency
histograms (FEHs) were employed as the equilibrium distributions of states.
Because of the uncertainty posed by experimental limitations (a relatively
small number of photons in time bins, and a long binning time in comparison to
the characteristic times of interconversion of protein molecules), it was not
quite clear whether these FEHs could be used as the PMFs. However, it has been
found that for the denaturant concentrations at which the FEH peak for the
target folded (unfolded) state is resolved, the rates constants of folding
(unfolding) change with the denaturant concentration in satisfactory agreement
with those extracted from photon trajectories with the use of theoretic models
gopich2009decoding ; gopich2010fret . Also, it has been found the FPT
distributions are single-exponential, which is consistent with the assumption
of two-state kinetics in the theoretical models gopich2009decoding ;
gopich2010fret  for extracting the rate constants from the FRET data. To
obtain absolute values of the rate constants, a time-dependent quantity is
required, which, on the one hand, would be known in the experiment (physical
or computational) and, on the other, could be obtained in the MC simulations
based on the given equilibrium distribution of states, e.g., the self-
diffusion coefficient, if this distribution were in physical space
huitema1999can . The promising feature of the present approach is that it is
based solely on the equilibrium distributions of states, without introducing
any additional parameters to perform simulations, which suggests its
applicability to other complex systems.

###### Acknowledgements.

I thank Irina Gopich for valuable discussions of the results of the FRET-
experiments for folding of $\alpha_{3}$D protein, and Hoi Sung Chung for
sharing the corresponding FRET-efficiency histograms. The work was supported
by a contract with the IT SB RAS.

## Author Declarations

Conflict of Interest
The authors have no conflicts to disclose.

## Data Availability Statement

AVAILABILITY OF DATA |  STATEMENT OF DATA AVAILABILITY
---|---
Data available in article or on request from the author |  The data that support the findings of this study are available within the article or from the author upon reasonable request.

## References

  * (1) B. Schuler, E. A. Lipman, and W. A. Eaton, Nature 419, 743 (2002).
  * (2) F. Huang, L. Ying, and A. R. Fersht, Proc. Natl. Acad. Sci. USA 106, 16239 (2009).
  * (3) H. S. Chung, I. V. Gopich, K. McHale, T. Cellmer, J. M. Louis, and W. A. Eaton, J. Phys. Chem. A 115, 3642 (2011).
  * (4) J. Liu, L. A. Camposa, M. Cerminara, X. Wang, R. Ramanathan, D. S. Englisha, and V. Muñoz, Proc. Natl. Acad. Sci. USA 109, 179 (2012).
  * (5) B. Schuler and H. Hofmann, Curr. Opin. Struct. Biol. 23, 36 (2013).
  * (6) W. A. Eaton, J. Phys. Chem. B 125, 3452 (2021).
  * (7) H. Yu, A. N. Gupta, X. Liu, K. Neupane, A. M. Brigley, I. Sosova, and M. T. Woodside, Proc. Natl. Acad. Sci. USA 109, 14452 (2012).
  * (8) M. Batchelor, K. Papachristos, M. Stofella, Z. T. Yew, and E. Paci, Biochim. Biophys. Acta 1864, 129613 (2020).
  * (9) C. Bustamante, L. Alexander, K. Maciuba, and C. M. Kaiser, Annu. Rev. Biochem. 89, 443 (2020).
  * (10) S. A. McKinney, C. Joo, and T. Ha, Biophys. J. 91, 1941 (2006).
  * (11) I. V. Gopich and A. Szabo, J. Phys. Chem. B 113, 10965 (2009).
  * (12) I. V. Gopich and A. Szabo, J. Phys. Chem. B 114, 15221 (2010).
  * (13) R. Ramanathan and V. Muñoz, J. Phys. Chem. B 119, 7944 (2015).
  * (14) I. S. Osad’ko, J. Phys. Chem. C 120, 697 (2016).
  * (15) G. Hummer and A. Szabo, Proc. Natl. Acad. Sci. USA 98, 3658 (2001).
  * (16) G. Hummer and A. Szabo, Acc. Chem. Res. 38, 504 (2005).
  * (17) C. Hyeon and D. Thirumalai, J. Phys.: Condens. Matter 19, 113101 (2007).
  * (18) O. K. Dudko, G. Hummer, and A. Szabo, Proc. Natl. Acad. Sci. USA 105, 15755 (2008).
  * (19) 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 330, 341 (2010).
  * (20) K. Lindorff-Larsen, S. Piana, R. O. Dror, and D. E. Shaw, Science 334, 517 (2011).
  * (21) C. Camilloni and F. Pietrucci, Advances in Physics: X 3, 1477531 (2018).
  * (22) Y. Sugita and Y. Okamoto, Chem. Phys. Lett. 314, 141 (1999).
  * (23) Y. M. Rhee and V. S. Pande, Biophys. J. 84, 775 (2003).
  * (24) F. Rao and A. Caflisch, J. Chem. Phys. 119, 4035 (2003).
  * (25) A. M. Ferrenberg and R. H. Swendsen , Phys. Rev. Lett. 63, 1195 (1989).
  * (26) J.-E. Shea and C. L. Brooks III, Annu. Rev. Phys. Chem. 52, 499 (2001).
  * (27) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. N. Teller, and E. Teller, J. Chem. Phys 21, 1087 (1953).
  * (28) D. Frenkel, Introduction to Monte Carlo Methods, Computational Soft Matter: From Synthetic Polymers to Proteins, Lecture Notes, John von Neumann Institute for Computing, Jülich, NIC Series, N. Attig, K. Binder, H. Grubmüller, K. Kremer (Eds.), 2004, 23, 29-60.
  * (29) S. F. Chekmarev, A. Yu. Palyanov, and M. Karplus, Phys. Rev. Lett. 100, 018107 (2008).
  * (30) I. V. Kalgin, S. F. Chekmarev, and M. Karplus, J. Phys. Chem. B 118, 4287 (2014).
  * (31) B. Roux, Comput. Phys. Comm. 91, 275 (1995).
  * (32) D. Trzesniak, A.-P. E. Kunz, and W. F. van Gunsteren, ChemPhysChem 8, 162 (2007).
  * (33) H. E. A. Huitema and J. P. van der Eerden, J. Chem. Phys. 110, 3267 (1999).
  * (34) G. Rutkai and T. Kristóf, J. Chem. Phys. 132, 104107 (2010).
  * (35) E. Sanz and D. Marenduzzo, J. Chem. Phys. 132, 194102 (2010).
  * (36) K. M. Bal and E. C. Neyts, J. Chem. Phys. 141, 204104 (2014).
  * (37) K. Kikuchi, M. Yoshida, T. Maekawa, and H. Watanabe, Chem. Phys. Lett. 185, 335 (1991).
  * (38) G. Tiana, L. Sutto, and R. A. Broglia, Physica A: Stat. Mech. Appl. 380, 241 (2007).
  * (39) N. Gō, Annu. Rev. Biophys. Bioeng. 12, 183 (1983).
  * (40) S. F. Chekmarev, J. Chem. Phys. 139, 145103 (2013).
  * (41) N. H. Andersen, K. A. Olsen, R. M. Fesinmeyer, X. Tan, F. M. Hudson, L. A. Eidenschink, and S. R. Farazi, J. Am. Chem. Soc. 128, 6101 (2006).
  * (42) T. X. Hoang and M. Cieplak, J. Chem. Phys. 112, 6851 (2000).
  * (43) R. Biswas and D. R. Hamann, Phys. Rev. B 34, 895 (1986).
  * (44) S. Miyazawa and R. L. Jernigan, J. Mol. Biol. 256, 623 (1996).
  * (45) S. F. Chekmarev, Phys. Chem. Chem. Phys. 23, 17856 (2021).
  * (46) R. B. Best, G. Hummer, and W. A. Eaton, Natl. Acad. Sci. USA 110, 17874 (2013).
  * (47) V, Muñoz, P. A. Thompson, J. Hofrichter, and W. A. Eaton, Nature 390, 196 (1997).
  * (48) R. Zhou, B. J. Berne, and R. Germain, Proc. Natl. Acad. Sci. USA 98, 14931 (2001).
  * (49) Y. Gao, Y. Li, L. Mou, B. Lin, J. Z. H. Zhang, and Y Mei, Sci. Rep. 5, 10359 (2015).
  * (50) A. R. Dinner, T. Lazaridis, and M. Karplus, Natl. Acad. Sci. USA 96, 9068 (1999).
  * (51) B. Zagrovic, E. J. Sorin, and V. Pande, J. Mol. Biol. 313, 151 (2001).