Entropy and stochastic properties in catalysis at nanoscale

The equations for the mean substrate concentration, $\psi ,$ and the mean concentration of enzyme-substrate complex, $\phi ,$ that describe the macroscopic part are:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d\psi }{dt}=-{k}_{1}^{{\prime} }\left({n}_{10}-\phi \right)\psi +{k}_{2}\phi \\ \displaystyle \frac{d\phi }{dt}={k}_{1}^{{\prime} }\left({n}_{10}-\phi \right)\psi -\left({k}_{2}+{k}_{3}\right)\phi \end{array}\end{eqnarray} \tag{ 17 }$$

This approach is versatile enough to be of use at any enzyme-substrate ratio, since these differential equations can be solved numerically for different initial conditions; but it will be used with the usual proportions of reactants found in in-vitro experiments, where the initial substrate concentration is much greater than the initial enzyme concentration. The logic behind this is that smaller quantities of enzyme can catalyze much larger quantities of substrate, thus it is a topic of efficiency; however, this is not the only situation that can be studied in the laboratories. Albe et al [30] report the progress of a chemical reaction with different enzyme-substrate ratios, including cases where these ratios are inverted, so that the amount of enzyme exceeds that of the substrate.

Equations (17) can take a very practical form when multiplied by $\tfrac{1}{{k}_{1}^{{\prime} }},$ and one can define the evolution parameter $s={k}_{1}^{{\prime} }t,$ with units of $\tfrac{1}{M}.$ The resulting expressions are:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d\psi }{ds}=-\left({n}_{10}-\phi \right)\psi +r\phi \\ \displaystyle \frac{d\phi }{dt}=\left({n}_{10}-\phi \right)\psi -{k}_{M}\phi \end{array}\end{eqnarray} \tag{ 18 }$$

Where $r={k}_{2}/{k}_{1}^{{\prime} }$ is an efficiency parameter that measures the rate of dissociation of the $ES$ complex that is not transformed into product.

The quasi-stationary state that is of interest in the in-vitro experiments is obtained through the condition that the density of enzyme-substrate complex changes very little across time, which is mathematically expressed as $\tfrac{d{\phi }_{s}}{dt}\cong 0.$ This is the topic of discussion in this subsection.

In terms of the densities, the conservation laws take the form of

$$\begin{eqnarray}&&\begin{array}{c}{n}_{10}=e\left(t\right)+\phi \left(t\right),\,{n}_{20}=\psi \left(t\right)+\phi \left(t\right)+P\left(t\right)\end{array}\end{eqnarray} \tag{ 19 }$$

Where $e\left(t\right)$ and $P\left(t\right)$ are the densities of the enzyme and the product, respectively. It is possible to demonstrate that the time evolution equations can be obtained through the following expressions:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{de}{dt}=-\displaystyle \frac{d\phi }{dt},\,\displaystyle \frac{dP}{dt}=-{k}_{3}\phi \left(t\right)\end{array}\end{eqnarray} \tag{ 20 }$$

Introducing the quasi-stationary condition in the time evolution equation of $\phi \left(t\right)$ in (17) results in

$$\begin{eqnarray*}&&{k}_{1}^{{\prime} }\left({n}_{10}-{\phi }_{s}\right){\psi }_{s}-\left({k}_{2}+{k}_{3}\right){\phi }_{s}\cong 0\end{eqnarray*}$$

The conservation laws yield the relation $\left({n}_{10}-{\phi }_{s}\right)={e}_{s},$ such that

$$\begin{eqnarray}&&\begin{array}{c}{\phi }_{s}\cong \displaystyle \frac{{k}_{1}^{{\prime} }{e}_{s}{\psi }_{s}}{{k}_{2}+{k}_{3}}=\displaystyle \frac{{e}_{s}{\psi }_{s}}{{k}_{M}}\end{array}\end{eqnarray} \tag{ 21 }$$

where ${k}_{M}=\tfrac{{k}_{2}+{k}_{3}}{{k}_{1}^{{\prime} }}$ and receives the well-known name of the Michaelis-Menten constant.

Denoting as $V$ the rate of growth of the product density $P\left(t\right),$ one obtains

$$\begin{eqnarray*}&&V=\displaystyle \frac{dP\left(t\right)}{dt}\end{eqnarray*}$$

The in-vitro experiments are performed with amounts on the order micromoles. In such cases the behavior observed is of a slow decrease of the functions $\psi \left(t\right)$ and $\phi \left(t\right).$ It is then where the condition of $\tfrac{d{\phi }_{s}}{dt}$ is applicable.

From the deterministic equation for $\tfrac{d\phi }{dt}$ in (18) results

$$\begin{eqnarray}&&\begin{array}{c}\left({n}_{10}-\phi \right)\psi -{k}_{M}\phi \cong 0\end{array}\end{eqnarray} \tag{ 22 }$$

One can obtain the following expression for the mean concentration of enzyme-substrate complex:

$$\begin{eqnarray*}&&{\phi }_{s}\cong \displaystyle \frac{{n}_{10}{\psi }_{s}}{{k}_{M}+{\psi }_{s}}\end{eqnarray*}$$

Using de definition $V=\tfrac{dP\left(t\right)}{dt}$ and the time evolution of $P\left(t\right)$ in (20), the demonstration of the following expression is straightforward.

$$\begin{eqnarray*}&&V\cong \displaystyle \frac{{V}_{\max }{\psi }_{s}}{{k}_{M}+{\psi }_{s}}\end{eqnarray*}$$

Where ${V}_{\max }={k}_{3}{n}_{10}.$ The usual notation is recovered by simply rewriting ${\psi }_{s}$ as $\left[S\right].$

It is common practice to indistinctively use the terms 'stationary' and 'equilibrium' as synonyms, but it is imperative to establish a clear distinction between them. This section is dedicated to this purpose, showing the difference between a stationary state achieved in in-vitro experiments and the equilibrium state. This latter one corresponds to the mathematical condition:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{d\psi }{dt}=0,\,\displaystyle \frac{d\phi }{dt}=0\end{array}\end{eqnarray} \tag{ 23 }$$

such that the following equations are followed:

$$\begin{eqnarray}&&\begin{array}{l}-{k}_{1}^{{\prime} }\left({n}_{10}-{\phi }_{e}\right){\psi }_{e}+{k}_{2}{\phi }_{e}=0\\ {k}_{1}^{{\prime} }\left({n}_{10}-{\phi }_{e}\right){\psi }_{e}-\left({k}_{2}+{k}_{3}\right){\phi }_{e}=0\end{array}\end{eqnarray} \tag{ 24 }$$

The algebraic solution is ${\psi }_{e}=0,\,{\phi }_{e}=0,$ it corresponds to the complete consumption of the substrate and enzyme-substrate complex. It is for this reason that the state of equilibrium is only attainable when all the substrate has been consumed, such that only the initial amounts of free enzymes and the product remain in the system.

It is possible to study the stability of the equilibrium state. Linearizing the differential equations around $\left({\psi }_{e},{\phi }_{e}\right)=\left(0,0\right)$ results in:

$$\begin{eqnarray*}&&\displaystyle \frac{d}{dt}\left(\begin{array}{c}\psi \left(t\right)\\ \phi \left(t\right)\end{array}\right)=\left(\begin{array}{ll}-{k}_{1}^{{\prime} }{n}_{10} & {k}_{2}\\ {k}_{1}^{{\prime} }{n}_{10} & -\left({k}_{2}+{k}_{3}\right)\end{array}\right)\left(\begin{array}{c}\psi \\ \phi \end{array}\right)\end{eqnarray*}$$

To ease reading, we introduce a shorter notation:

$$\begin{eqnarray*}&&A={k}_{1}^{{\prime} }{n}_{10}\,,\,C={k}_{2},\,D={k}_{3},\,R=\sqrt{-4AD+{\left(A+C+D\right)}^{2}}\end{eqnarray*}$$

Thus, the eigenvalues of the matrix above are given as:

$$\begin{eqnarray*}&&{\lambda }_{1,2}=-\displaystyle \frac{1}{2}\left(A+C+D\pm R\right)\end{eqnarray*}$$

The transition rates and initial conditions are defined positive, this results in the eigenvalues satisfying the inequality: ${\lambda }_{1,2}\lt 0.$ The fundamental solutions $\left\{{e}^{{\lambda }_{1}t},{e}^{{\lambda }_{2}t}\right\}$ generate the most general solution:

$$\begin{eqnarray*}&&\left(\begin{array}{l}\psi \left(t\right)\\ \phi \left(t\right)\end{array}\right)=\left(\begin{array}{l}A{e}^{{\lambda }_{1}t}+B{e}^{{\lambda }_{2}t}\\ E{e}^{{\lambda }_{1}t}+F{e}^{{\lambda }_{2}t}\end{array}\right)\end{eqnarray*}$$

Which will always converge towards the equilibrium state. This mathematic expression explains the profile of the curve for the $ES$ complex seen in figure 1 and the decreasing stage shown in figure 14.

The following equation is called the Fokker-Planck equation (FPE):

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\partial P\left(\vec{q},t\right)}{\partial t}=\displaystyle \frac{\partial \left[{A}_{\mu }\left(\vec{q},\psi ,\phi \right)P\left(\vec{q},t\right)\right]}{\partial {q}_{\mu }}+\displaystyle \frac{1}{2}\displaystyle \sum _{\mu =1}^{2}\displaystyle \sum _{\nu =1}^{2}{D}_{\mu \nu }\left(\psi ,\phi \right)\displaystyle \frac{{\partial }^{2}P\left(\vec{q},t\right)}{\partial {q}_{\mu }{q}_{\nu }}\end{array}\end{eqnarray} \tag{ 25 }$$

The expression above describes the probability that the fluctuation of substrate concentration takes the value ${q}_{1},$ and the fluctuation of enzyme-substrate complex concentration takes the value ${q}_{2},$ at time $t.$ The systematic term can be written as

$$\begin{eqnarray}&&\begin{array}{c}{A}_{\mu }\left(\vec{q},\psi ,\phi \right)=\displaystyle \sum _{\nu =1}^{2}{L}_{\mu \nu }\left(\psi ,\phi \right){q}_{\nu }\end{array}\end{eqnarray} \tag{ 26 }$$

Where ${A}_{\mu }$ is the flux of random concentration fluctuations of $S$ and $ES.$

And its general solution [31] is a gaussian function of the form:

$$\begin{eqnarray}&&\begin{array}{c}P\left(\vec{q},t\right)=\displaystyle \frac{1}{2\pi \sqrt{Det\left({\rm{\Xi }}\left(t\right)\right)}}{e}^{-\displaystyle \frac{1}{2}\left({q}_{\mu }-\left\langle {q}_{\mu }\left(t\right)\right\rangle \right){{\rm{\Xi }}}^{-1}\left(t\right)\left({q}_{\nu }-\left\langle {q}_{\nu }\left(t\right)\right\rangle \right)}\end{array}\end{eqnarray} \tag{ 27 }$$

with ${\rm{\Xi }}\left(t\right)=\left\langle {q}_{\mu }{q}_{\nu }\right\rangle -\left\langle {q}_{\mu }\right\rangle \left\langle {q}_{\nu }\right\rangle$ as the self-correlation matrix, and ${{\rm{\Xi }}}^{-1}\left(t\right)$ as its inverse. This is the matrix that contains the variance of the random concentration fluctuations of $E$ and $ES,$ as well as their correlations.

In the case under study, the elements of $\left\{{L}_{\mu \nu }\right\}$ and $\left\{{D}_{\mu \nu }\right\}$ are given by:

$$\begin{eqnarray}&&\begin{array}{l}{L}_{11}=-{k}_{1}^{{\prime} }\left({n}_{10}-\phi \right),\,{L}_{12}={k}_{1}^{{\prime} }\psi +{k}_{2}\\ {L}_{21}={k}_{1}^{{\prime} }\left({n}_{10}-\phi \right)\,,\,{L}_{22}=-{k}_{1}^{{\prime} }\psi -{k}_{2}-{k}_{3}\end{array}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\begin{array}{l}{D}_{11}={k}_{1}^{{\prime} }\left({n}_{10}-\phi \right)\psi +\left(1+{k}_{2}\right)\phi \\ {D}_{22}={k}_{1}^{{\prime} }\left({n}_{10}-\phi \right)\psi +\left({k}_{2}+{k}_{3}\right)\phi \\ {D}_{12}=-{k}_{1}^{{\prime} }\left({n}_{10}-\phi \right)\psi -{k}_{2}\phi \end{array}\end{eqnarray} \tag{ 29 }$$

Given a physical magnitude $g(\vec{q},t)$ that is relevant to the system, its mean can be calculated by:

$$\begin{eqnarray}&&\begin{array}{c}\left\langle g\left(\vec{q},t\right)\right\rangle \left(t\right)=\displaystyle {\int }_{-\infty }^{\infty }\displaystyle {\int }_{-\infty }^{\infty }g\left(\vec{q},t\right)P\left(\vec{q},t\right)d{q}_{1}d{q}_{2}\end{array}\end{eqnarray} \tag{ 30 }$$

The time evolution equations of the mean of the fluctuations of the concentrations of $E$ and $ES$ take the form of:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{d}{dt}\left(\begin{array}{c}\left\langle {q}_{1}\right\rangle \left(t\right)\\ \left\langle {q}_{2}\right\rangle \left(t\right)\end{array}\right)=\left(\begin{array}{cc}{L}_{11} & {L}_{12}\\ {L}_{21} & {L}_{22}\end{array}\right)\left(\begin{array}{c}\left\langle {q}_{1}\right\rangle \left(t\right)\\ \left\langle {q}_{2}\right\rangle \left(t\right)\end{array}\right)\end{array}\end{eqnarray} \tag{ 31 }$$

and the time evolution equations for the self-correlation functions are:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{d}{dt}\left(\begin{array}{c}{{\rm{\Xi }}}_{11}\left(t\right)\\ {{\rm{\Xi }}}_{12}\left(t\right)\\ {{\rm{\Xi }}}_{22}\left(t\right)\end{array}\right)=\left(\begin{array}{lll}2{L}_{11} & 2{L}_{12} & 0\\ {L}_{21} & {L}_{11}+{L}_{22} & {L}_{12}\\ 0 & 2{L}_{21} & 2{L}_{22}\end{array}\right)\left(\begin{array}{c}{{\rm{\Xi }}}_{11}\left(t\right)\\ {{\rm{\Xi }}}_{12}\left(t\right)\\ {{\rm{\Xi }}}_{22}\left(t\right)\end{array}\right)+\left(\begin{array}{c}{D}_{11}\left(t\right)\\ {D}_{12}\left(t\right)\\ {D}_{22}\left(t\right)\end{array}\right)\end{array}\end{eqnarray} \tag{ 32 }$$

Defining $\vec{{\rm{\Xi }}}\left(t\right)=\left[{{\rm{\Xi }}}_{11}\left(t\right),{{\rm{\Xi }}}_{12}\left(t\right),{{\rm{\Xi }}}_{22}\left(t\right)\right]$ and $\vec{D}\left(t\right)=\left[{D}_{11}\left(t\right),{D}_{12}\left(t\right),{D}_{22}\left(t\right)\right],$ the equation above can be rewritten as follows:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{d\vec{{\rm{\Xi }}}\left(t\right)}{dt}=M\left(t\right)\vec{{\rm{\Xi }}}\left(t\right)+\vec{D}\left(t\right)\end{array}\end{eqnarray} \tag{ 33 }$$

where

$$\begin{eqnarray}\begin{array}{c}M\left(t\right)=\left(\begin{array}{lll}2{L}_{11}\left(t\right) & 2{L}_{12}\left(t\right) & 0\\ {L}_{12}\left(t\right) & {L}_{11}\left(t\right)+{L}_{22}\left(t\right) & {L}_{12}\left(t\right)\\ 0 & 2{L}_{21}\left(t\right) & 2{L}_{22}\left(t\right)\end{array}\right)\end{array}\end{eqnarray} \tag{ 34 }$$

Defining the entropy of fluctuation of the concentrations of $E$ and $ES$ as

$$\begin{eqnarray}&&\begin{array}{c}S\left(t\right)=-\left\langle \mathrm{ln}\,P\left(\vec{q},t\right)\right\rangle =-\displaystyle {\int }_{-\infty }^{\infty }\displaystyle {\int }_{-\infty }^{\infty }P\left(\vec{q},t\right)\,\mathrm{ln}\left[P\left(\vec{q},t\right)\right]d{q}_{1}d{q}_{2}\end{array}\end{eqnarray} \tag{ 35 }$$

Substituting the general solution in the expression for the entropy (35), results in expression (36) seen below:

$$\begin{eqnarray}&&\begin{array}{c}S\left(t\right)=\displaystyle \frac{1}{2}\,\mathrm{ln}\left[Det\left({\rm{\Xi }}\left(t\right)\right)\right]+\,\mathrm{ln}\left(2\pi e\right)\end{array}\end{eqnarray} \tag{ 36 }$$

Given two temporal points such that ${t}_{1}\lt {t}_{2},$ the change in entropy is

$$\begin{eqnarray}&&\begin{array}{c}{\rm{\Delta }}S=S\left({t}_{2}\right)-S\left({t}_{1}\right)=\displaystyle \frac{1}{2}\,\mathrm{ln}\left[\displaystyle \frac{Det\left({\rm{\Xi }}\left({t}_{2}\right)\right)}{Det\left({\rm{\Xi }}\left({t}_{1}\right)\right)}\right]\end{array}\end{eqnarray} \tag{ 37 }$$

so, the condition of ${\rm{\Delta }}S\geqslant 0$ translates into the following condition:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{Det\left({\rm{\Xi }}\left({t}_{2}\right)\right)}{Det\left({\rm{\Xi }}\left({t}_{1}\right)\right)}\gt 1\end{array}\end{eqnarray} \tag{ 38 }$$

Expression (38) will later be utilized to confirm that the entropy of fluctuation decreases when the reaction is taking place. Expression (27) provides the mathematical form of the time dependent probability density for the fluctuation variables which, as will be shown in the numerical example in a later section, performs a clockwise turning motion as the reaction progresses. This result showcases the importance of fluctuations in the dynamics of catalysis; therefore, in-depth study of this topic is relevant.

If the phenomenon is analyzed using the dynamics of stochastic velocities, the behavior observed can be better understood. Thus, the section below is dedicated to the introduction of this formalism that may prove useful to the uninitiated reader.

The time evolution of the fluctuations can be described by the motion of a state point $\vec{q}.$ It is common to assume that its study is completed once its behavior has been formulated, as we did in the previous section. We will see that it is possible to add new knowledge to better understand the behavior of the state point $\vec{q}.$

In this section we present an intuitive approach to stochastic velocities. We begin with an analysis of the difficulty found when attempting to define velocities in stochastic processes. The typical example given of a stochastic process is the Brownian motion; this describes the conduct of a micrometric mass floating inside a liquid at temperature $T.$ Seen through a modern video camera, its movement takes place in two dimensions, but for ease of writing mathematical expressions, we consider only one dimension. The model treats the Brownian particle as if it is a point mass, and due to the movement possessing random behavior, each point $z$ has a time dependent probability associated to it. This probability is a probability density function, denoted as $\rho \left(z,t\right),$ such that a given interval, $\left(z,z+dz\right),$ on the line of accessible positions, the expression $\rho \left(z,t\right)dz$ yields the probability of the particle being found within that interval at time $t.$

It can be demonstrated that $\rho \left(z,t\right)$ follows an equation of the form shown below:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\partial \rho \left(z,t\right)}{\partial t}=D\displaystyle \frac{{\partial }^{2}\rho \left(z,t\right)}{\partial {z}^{2}}\end{array}\end{eqnarray} \tag{ 39 }$$

where $D$ is named the diffusion coefficient. With the initial condition $\rho \left(z,t=0\right)=\delta \left(z-{z}_{0}\right)$ the solution obtained for (39) is

$$\begin{eqnarray}&&\begin{array}{c}\rho \left(z,t\right)=\displaystyle \frac{1}{\sqrt{4\pi Dt}}{e}^{-\displaystyle \frac{{\left(z-{z}_{0}\right)}^{2}}{4\pi Dt}}\end{array}\end{eqnarray} \tag{ 40 }$$

In this case, the probability densities evolve with time. From (40), the statistical properties $\left\langle z\right\rangle =0$ and $\left\langle {z}^{2}\right\rangle =2Dt$ can be demonstrated. It is important to note that the mean movement is zero and that the standard deviation changes with time as ${t}^{1/2}.$ In the approach presented by Paul Langevin in 1908, the Second Law of Newton is applied to a Brownian particle of mass $m$ with a force $F\left(x\right)$ acting upon it; plus a friction force proportional to the velocity, $-\beta v;$ plus a stochastic force, $\eta \left(t\right)$ with the following properties:

$$\begin{eqnarray}&&\begin{array}{c}\left\langle \eta \left(t\right)\right\rangle =0,\,\left\langle \eta \left(t\right)\eta \left(s\right)\right\rangle =C\delta \left(t-s\right)\end{array}\end{eqnarray} \tag{ 41 }$$

where $C$ is a constant. There are technical reasons for calling white noise a random magnitude that displays these properties. The fundamental dynamic law resulting from this is called the Langevin equation, and is written as shown below:

$$\begin{eqnarray}&&\begin{array}{c}m\displaystyle \frac{dv}{dt}=F\left(x\right)-\beta v+\eta \left(t\right)\end{array}\end{eqnarray} \tag{ 42 }$$

Using the method of moments, one obtains similar results for the Brownian motion. If the trajectories as functions of time $t$ are to be plotted, the resulting graphs would present functions with very sharp peaks and valleys, therefore one can perceive intuitively that there would be problems when attempting to define the displacement velocity as $v=\tfrac{dx}{dt}.$ Next, we will focus on this dilemma with close attention.

In mathematical terms, the Weiner process has been defined to formalize the events that follow the Brownian motion. Taking as a starting point the case where $D=1$ and denoting the process as $W\left(t\right),$ the postulated properties are:

• 1.  
  $W\left(t=0\right)=0.$
• 2.  
  $W\left(t\right)$ is continuous.
• 3.  
  $W\left(t\right)$ has independent increments, that is to say, if it follows that $0\leqslant {s}_{1}\leqslant {t}_{1}\leqslant {s}_{2}\leqslant {t}_{2},$ then the differences $W\left({t}_{1}\right)=W\left({s}_{1}\right)$ and $W\left({t}_{2}\right)-W\left({s}_{2}\right)$ are independent random variables.
• 4.  
  $W\left(t-s\right)$ is a random variable with normal distribution of mean $\mu =0$ and variance ${\sigma }^{2}=t-s.$

In rigorous terms, equation (42) presents a difficulty that is analyzed next. Using the traditional concepts of defining the rate of change in a trajectory $z\left(t\right),$ we take the quotient of finite increments:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{{\rm{\Delta }}z}{{\rm{\Delta }}t}=\displaystyle \frac{z\left(t+h\right)-z\left(t\right)}{h}\end{array}\end{eqnarray} \tag{ 43 }$$

but from the fourth postulate one must have $z\left(t+h\right)-z\left(t\right)=\sqrt{h},$ where it follows that $\tfrac{{\rm{\Delta }}z}{{\rm{\Delta }}t}=\tfrac{\sqrt{h}}{h}.$ Therefore, the quotient diverges in the limit $h\to 0;$ consequently, it is impossible to define the rate of change through traditional methods.

For the issue encountered above, instead of using the Langevin equation, as it was originally formulated, it is better to consider an approach using finite differences to suggest an equation that avoids the use of derivatives and translates all calculations to integrals, giving place to two types of integral calculus: that of Kiyoshi Ito and of Ruslan Stratonovich [32]. A less known option is the one developed by Edward Nelson, based in a system of averages over the realizations of the stochastic processes, and gives place to the concept of stochastic velocities [22, 33]. These have been used to describe quantum phenomena from a stochastic perspective, giving rise to a line of work called stochastic mechanics. In this work we take advantage of the mathematical tools developed and use them in our topic of interest. The condition is that the stochastic process must be describable by means of Fokker-Planck equations (FPE) [25].

When an ensemble is associated to a stochastic process $\vec{x}=\left({x}_{1},\ldots ,{x}_{n}\right),$ a state space of dimension $n$ is available in which a point at a time $t$ corresponds to each member of the ensemble. A large number of members of the ensemble produce a cloud of points that move at random as time progresses; this idea allows the introduction of an analogy with a gas, such that it is possible to include in the description some properties typical of fluids. One of these is that of vorticity, which lets us know if the cloud of state points tends to rotate, or if it behaves like a fluid whose velocity is irrotational. This gives us the opportunity of using this concept to analyze the manner with which the agitation of this cloud of state points occurs, and with it establish a difference between a stationary process and that of a process at equilibrium. The former abides to the condition that the statistical moment of order $m$ must be time invariant:

$$\begin{eqnarray}&&\begin{array}{c}\left\langle {x}_{1}\left(t\right){x}_{2}\left(t\right)\ldots {x}_{m}\left(t\right)\right\rangle =\left\langle {x}_{1}\left(t+\tau \right){x}_{2}\left(t+\tau \right)\ldots {x}_{3}\left(t+\tau \right)\right\rangle \end{array}\end{eqnarray} \tag{ 44 }$$

while the latter also follows the condition that each of the possible transitions must be balanced out by a transition that occurs on the opposite side. This condition is given the name of detailed balance. If the detailed balance is not present, then the cloud of state points tends to rotate, therefore, the use of the concept of vorticity contributes to the understanding of the dynamics of the gas cloud. Vorticity is defined as the curl of the velocity thus it is necessary to revise this point.

Smoothening the trajectory using a moving average

The mean over the realizations, as originally developed by Nelson, can be understood with the concept of moving averages used in the study of time series. In this section we introduce these concepts.

Let $\vec{x}\left(t\right)$ be a stochastic process that occurs inside a state space ${\mathscr{U}},$ and let the points $\left\{\vec{x}^{\prime} ,\vec{x},\vec{x}^{\prime\prime} \right\}\in {\mathscr{U}},$ such that they are reached by some realizations of the stochastic process $\vec{x}$ at times ${t}^{{\prime} },t,t^{\prime\prime} ,$ respectively, where $t^{\prime} \lt t\lt t^{\prime\prime} .$ Let the increments be ${{\rm{\Delta }}}_{+}\vec{x}=\vec{x}^{\prime\prime} -\vec{x},\,{{\rm{\Delta }}}_{-}\vec{x}=\vec{x}-\vec{x}^{\prime} ,$ see figure 5.

Zoom In Zoom Out Reset image size

Next, we will explain the relations that must be followed between the time intervals involved in the description; for that purpose, we continue using the Brownian movement as a case study. Suppose a camera capable of registering the random movements of a state point; due to the difference in mass, the collision of a single molecule against the Brownian particle does not produce an effect that is registerable by the lab instrument. What moves the Brownian particle is the difference in the number of collisions that it receives because of the random fluctuations of the density of particles comprising the surrounding medium. By this manner, the path traced by the Brownian particle, as observed by a camera recording through an optical microscope, are polygonal shaped. However, because of the easiness in the mathematical language, we say that a collision occurs each time there is a change in the path taken by the particle of interest.

There are three instants of time that are relevant in this approach:

• 1.  
  ${t}_{c},$ the time required to accumulate enough collisions capable of causing a registerable random change in the path of the particle.
• 2.  
  ${T}_{0},$ the time that passes between two successive frames captured by the camera. If ${T}_{0}$ is too short, the displacements registered may follow the relation ${\rm{\Delta }}x\sim t,$ as is the case in classical mechanics; but it is more common to find $\left\langle {\rm{\Delta }}x\right\rangle \sim {t}^{1/2}$ in the Brownian case. For the latter conduct it is necessary that ${T}_{0}\gg {t}_{c}.$
• 3.  
  ${\rm{\Delta }}$t, the time necessary to smoothen trajectories by the moving averages method. To study the movement of the state points, the stochastic trajectories are smoothed by introducing moving averages; these are calculated in time intervals ${\rm{\Delta }}t$ that must be long enough to include various spikes of the trajectory, as can be seen in figure 7, but also short enough that the camera registering the data cannot tell that the trajectory has been smoothed.

Therefore, ${\rm{\Delta }}t\gg {T}_{0}\gg {t}_{c}$ must be followed. This regime is called the coarse grain time scale, as shown in figure 6 below.

Zoom In Zoom Out Reset image size

The moving average of $p$ points is defined as ${x}_{k}=\tfrac{1}{p}\displaystyle {\sum }_{i=k-\tfrac{p}{2}}^{k+\tfrac{p}{2}}{x}_{i},$ where $i$ takes values such that the sum is not out of the bounds of the interval. For continuous signals, the moving average is defined as: ${\left\langle x\right\rangle }_{{\rm{\Delta }}t}\left(t\right)\equiv \displaystyle \frac{1}{{\rm{\Delta }}t}\displaystyle {\int }_{t^{\prime} }^{t^{\prime} +{\rm{\Delta }}t}x\left(s\right)ds.$ We now set out to study the derivative of a function $g\left(\vec{x},t\right).$ For this purpose, we establish a set of times expressed in (45) below:

$$\begin{eqnarray}&&\begin{array}{c}{t}_{1}\lt {t}_{2}\lt \ldots \lt {t}_{M}\end{array}\end{eqnarray} \tag{ 45 }$$

and considering finite increments of the function $g(\vec{x},t):$ ${\rm{\Delta }}{g}_{k}=g\left(\vec{x},{t}_{k}+1\right)-g(\vec{x},{t}_{k}),$ with $k=1,\ldots ,M-1,$ one can write:

$$\begin{eqnarray}&&\begin{array}{c}g\left(\vec{x},{t}_{M}\right)=g\left(\vec{x},{t}_{1}\right)+{\rm{\Delta }}{g}_{1}+{\rm{\Delta }}{g}_{2}+\ldots +{\rm{\Delta }}{g}_{M-1}\end{array}\end{eqnarray} \tag{ 46 }$$

Rearranging and multiplying by $\tfrac{1}{2{\rm{\Delta }}t}$ gives

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{g\left(\vec{x},{t}_{M}\right)-g\left(\vec{x},{t}_{1}\right)}{2{\rm{\Delta }}t}=\displaystyle \frac{1}{2{\rm{\Delta }}t}\displaystyle {\int }_{{t}^{^{\prime} }}^{t^{\prime} +2{\rm{\Delta }}t}\displaystyle \frac{\partial g\left(\vec{x},s\right)}{\partial s}\end{array}\end{eqnarray} \tag{ 47 }$$

This expression receives the name of coarse grain time derivative. Notice the incorporation of the sum of finite differences inside a moving average within an interval of width $2{\rm{\Delta }}t.$ Once the smoothening process has been applied, the resulting trajectories can be studied using the usual tools of calculus. An illustrative example is given in figure 7.

Zoom In Zoom Out Reset image size

The study of smoothened stochastic functions

Points in the state space and a statistical description of their movement

Suppose a statistical ensemble of equally prepared experiments at a macro scale. Each of these carry out a realization of the stochastic process $\vec{x}.$ At a given point in time, we have a cloud of points whose number is theoretically infinite. If one is to let the clock run out, a static image (like that of a photograph) would be substituted for a series of images of points moving at random, similar to a swarm of mosquitos swirling in summer. These will enter and exit of the previously marked region ${\mathcal B} ,$ as seen in figure 5. The question that follows is: How many points enter or exit ${\mathcal B}$ in a second? This problem is similar as counting the number of smoke particles in a given region of space.

The concept of systematic velocity, which has already been presented, measure the net number of state points that cross ${\mathcal B}$ per second.

Systematic derivative and systematic velocity

We define the systematic derivative as the mean over the ensemble of all the realizations

$$\begin{eqnarray}&&\begin{array}{c}{D}_{c}g\left(\vec{x},t\right)\equiv \displaystyle \frac{\left\langle g\left(\vec{x},t^{\prime\prime} \right)-g\left(\vec{x},t^{\prime} \right)\right\rangle }{2{\rm{\Delta }}t}\end{array}\end{eqnarray} \tag{ 48 }$$

To have an analytical representation of the previous definition, we introduce various hypotheses that lead to the Taylor series expansion. The hypotheses are:

• The stochasticity of the physical phenomenon is introduced by a source that is: stationary, isotropic and homogenous.
• The second statistical moments of the increments ${{\rm{\Delta }}}_{+}\vec{x}$ and ${{\rm{\Delta }}}_{-}\vec{x}$ are such that they follow:

$$\begin{eqnarray}&&\begin{array}{l}\left\langle {{\rm{\Delta }}}_{+}{x}_{i}{{\rm{\Delta }}}_{+}{x}_{j}\right\rangle =\left\langle {{\rm{\Delta }}}_{-}{x}_{i}{{\rm{\Delta }}}_{-}{x}_{j}\right\rangle \\ \left\langle {{\rm{\Delta }}}_{+}{x}_{i}{{\rm{\Delta }}}_{-}{x}_{j}\right\rangle =0\end{array}\end{eqnarray} \tag{ 49 }$$

The diffusion matrix can be defined as

$$\begin{eqnarray}&&\begin{array}{c}{D}_{ij}\equiv \displaystyle \frac{1}{2{\rm{\Delta }}t}\left\langle {{\rm{\Delta }}}_{+}x{{\rm{\Delta }}}_{+}{x}_{j}\right\rangle \end{array}\end{eqnarray} \tag{ 50 }$$

or also as

$$\begin{eqnarray}&&\begin{array}{c}{D}_{ij}\equiv \displaystyle \frac{1}{2{\rm{\Delta }}t}\left\langle {{\rm{\Delta }}}_{-}{x}_{i}{{\rm{\Delta }}}_{-}{x}_{j}\right\rangle \end{array}\end{eqnarray} \tag{ 51 }$$

Notice that the increments grow as ${\left({\rm{\Delta }}t\right)}^{1/2}$ so that the average that appears in the definition above grows as $\left({\rm{\Delta }}t\right).$ It also must be noted that the products of increments appear as coefficients in the second order derivatives of any function $d\left(x,t\right)$ expanded using Taylor series.

To work at higher orders of (${\rm{\Delta }}t$) would involve considering Taylor expansions with derivatives of orders of $3,4,\ldots ;$ although it is mathematically possible, there exist a restriction when the problem is translated to determining the probability $P(\vec{x},t)$ by means of a partial differential equation. The function $P(\vec{x},t)$ that satisfies an equation of order higher than 2 ceases to be defined as nonnegative for all $\vec{x}$ and $t,$ therefore it cannot be interpreted as a probability density function.

In vector calculus, a function $f\left(\vec{x},t\right)$ has a total time derivative that is of the form:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{df\left(\vec{x},t\right)}{dt}=\displaystyle \frac{\partial f\left(\vec{x},t\right)}{\partial t}+\displaystyle \sum _{i=1}^{n}\displaystyle \frac{\partial f}{\partial {x}_{i}}\displaystyle \frac{d{x}_{i}}{dt}\\ \displaystyle \frac{df\left(\vec{x},t\right)}{dt}=\displaystyle \frac{\partial f(\vec{x},t)}{\partial t}+{\vec{v}}_{c}\cdot {\rm{\nabla }}f\left(\vec{x},t\right)\end{array}\end{eqnarray} \tag{ 52 }$$

where ${\vec{v}}_{c}=\tfrac{d\vec{x}}{dt}$ is the velocity. The expression above is called the convective derivative.

This concept can be adapted to the case where the function $g\left(\vec{x},t\right)$ depends on a stochastic process $\vec{x}.$ This is the systematic derivative:

$$\begin{eqnarray}&&\begin{array}{c}{D}_{c}g\left(\vec{x},t\right)=\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}+{\vec{v}}_{c}\cdot {\rm{\nabla }}g\left(\vec{x},t\right)\end{array}\end{eqnarray} \tag{ 53 }$$

such that

$$\begin{eqnarray}&&\begin{array}{c}{\vec{v}}_{c}={D}_{c}\vec{x}\left(t\right)\equiv \displaystyle \frac{\left\langle \vec{x}\left({t}^{{\prime \prime} }\right)-\vec{x}\left({t}^{^{\prime} }\right)\right\rangle }{2{\rm{\Delta }}t}\end{array}\end{eqnarray} \tag{ 54 }$$

Access velocities, of exit and of diffusion

In fluid dynamics appears the phenomenon of swirls, or eddies, that cannot be understood with the systematic velocity alone; an illustrative example would be the that of cigar smoke traveling upwards through the air. To approach this topic let us consider a point $\vec{x}$ in the state space and both displacements, ${{\rm{\Delta }}}_{+}\vec{x}$ and ${{\rm{\Delta }}}_{-}\vec{x},$ as shown in figure 5. We have the next relations:

$$\begin{eqnarray}&&\begin{array}{c}\vec{x}^{\prime\prime} =\vec{x}+{{\rm{\Delta }}}_{+}\vec{x}\,,\,\vec{x}^{\prime} =\vec{x}+{{\rm{\Delta }}}_{-}\vec{x}\end{array}\end{eqnarray} \tag{ 55 }$$

If $\vec{x}$ is a state point contained in vicinity ${\mathcal B} ,$ it is understood that the displacement towards $\vec{x}^{\prime\prime} ,$ in a time interval ${\rm{\Delta }}t,$ is related to the exit of state points. Likewise, the displacement from $\vec{x}^{\prime}$ towards $\vec{x},$ also in a time interval ${\rm{\Delta }}t,$ is related with the entry of state points into vicinity ${\mathcal B} .$

Let us suppose we track a state point whose route is as follows:

• At instant $t-{\rm{\Delta }}t$ the state point is located at $\vec{x}^{\prime} .$
• At instant $t$ the state point is located at $\vec{x}^{\prime\prime} .$

Separating by components and working them individually, one can obtain

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{+}{x}_{i}={x}_{i}^{{\prime\prime} }-{x}_{i}\,{{\rm{\Delta }}}_{-}{x}_{i}={x}_{i}^{{\prime} }-{x}_{i}\end{eqnarray} \tag{ 56 }$$

We now consider the physical magnitude of the system denoted as $g(\vec{x}^{\prime\prime} ).$ To treat the displacement from $\vec{x}$ to $\vec{x}^{\prime\prime} ,$ the function is expanded in Taylor series as shown below:

$$\begin{eqnarray}&&\begin{array}{l}g\left(\vec{x}^{\prime\prime} ,t+{\rm{\Delta }}t\right)=g\left(\vec{x}+{{\rm{\Delta }}}_{+}\vec{x},t+{\rm{\Delta }}t\right)\\ \,=g\left(\vec{x},t\right)+\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}{\rm{\Delta }}t+\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}\left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\\ \,+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\left({x}_{j}^{{\prime\prime} }-{x}_{j}\right)+\ldots \end{array}\end{eqnarray} \tag{ 57 }$$

such that repeated indexes indicate a sum.

In the same fashion, one can study the displacement from $\vec{x}^{\prime}$ to $\vec{x},$ resulting in:

$$\begin{eqnarray}&&\begin{array}{l}g\left(\vec{x}^{\prime} ,t-{\rm{\Delta }}t\right)=g\left(\vec{x}+{{\rm{\Delta }}}_{-}\vec{x},t-{\rm{\Delta }}t\right)=g\left(\vec{x},t\right)-\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}{\rm{\Delta }}t\\ \,-\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}\left({x}_{i}-{x}_{i}^{{\prime} }\right)+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\left({x}_{i}-{x}_{i}^{{\prime} }\right)\left({x}_{j}-{x}_{j}^{{\prime} }\right)+\ldots \end{array}\end{eqnarray} \tag{ 58 }$$

The difference between $g\left(\vec{x}^{\prime\prime} \right)$ and $g\left(\vec{x}^{\prime} \right)$ is:

$$\begin{eqnarray}&&\begin{array}{l}g\left(\vec{x}^{\prime\prime} \right)-g\left(\vec{x}^{\prime} \right)=\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}2{\rm{\Delta }}t+\displaystyle \frac{\partial g}{\partial {x}_{i}}\left[{x}_{i}^{{\prime\prime} }-{x}_{i}+{x}_{i}-{x}_{i}^{{\prime} }\right]\\ \,+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}g}{\partial {x}_{i}\partial {x}_{j}}\left[\left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\left({x}_{j}^{{\prime\prime} }-{x}_{j}\right)-\left({x}_{i}-{x}_{i}^{{\prime} }\right)\left({x}_{j}-{x}_{j}^{{\prime} }\right)\right]+\ldots \end{array}\end{eqnarray} \tag{ 59 }$$

It is possible to find that the next relations are followed:

$$\begin{eqnarray}&&\begin{array}{l}{x}_{i}^{{\prime\prime} }-{x}_{i}+{x}_{i}-{x}_{i}^{{\prime} }={{\rm{\Delta }}}_{+}{x}_{i}+{{\rm{\Delta }}}_{-}{x}_{i}\\ \left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\left({x}_{j}^{{\prime\prime} }-{x}_{j}\right)-\left({x}_{i}-{x}_{i}^{{\prime} }\right)\left({x}_{j}-{x}_{j}^{{\prime} }\right)={{\rm{\Delta }}}_{+}{x}_{i}{{\rm{\Delta }}}_{+}{x}_{j}-{{\rm{\Delta }}}_{-}{x}_{i}{{\rm{\Delta }}}_{-}{x}_{j}\end{array}\end{eqnarray} \tag{ 60 }$$

Such that (59) can be written as:

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}g\left(\vec{x}^{\prime\prime} \right)-g\left(\vec{x}^{\prime} \right)=\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}2{\rm{\Delta }}t+\displaystyle \frac{\partial g}{\partial {x}_{i}}\left[{{\rm{\Delta }}}_{+}{x}_{i}+{{\rm{\Delta }}}_{-}{x}_{i}\right]\\ \,+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\left[{{\rm{\Delta }}}_{+}{x}_{i}{{\rm{\Delta }}}_{+}{x}_{j}-{{\rm{\Delta }}}_{-}{x}_{i}{{\rm{\Delta }}}_{-}{x}_{j}\right]+\ldots \end{array}\end{array}\end{eqnarray} \tag{ 61 }$$

Multiplying by $\tfrac{1}{2{\rm{\Delta }}t}$ and calculating the mean results in

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}\displaystyle \frac{\left\langle g\left(\vec{x}^{\prime\prime} \right)-g\left(\vec{x}^{\prime} \right)\right\rangle }{2{\rm{\Delta }}t}=\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}+\displaystyle \frac{\partial g}{\partial {x}_{i}}\displaystyle \frac{\left\langle {x}_{i}^{{\prime\prime} }-{x}_{i}^{{\prime} }\right\rangle }{2{\rm{\Delta }}t}\\ \,+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\displaystyle \frac{\left\langle {{\rm{\Delta }}}_{+}{x}_{i}{{\rm{\Delta }}}_{+}{x}_{j}-{{\rm{\Delta }}}_{-}{x}_{i}{{\rm{\Delta }}}_{-}{x}_{j}\right\rangle }{2{\rm{\Delta }}t}+\ldots \end{array}\end{array}\end{eqnarray} \tag{ 62 }$$

From (64) we have that the coefficient of the first derivative with respect to the position, is the ith component of the systematic velocity, ${v}_{i},$ which has been previously found. We also find the following:

$$\begin{eqnarray}&&\begin{array}{c}\left\langle {{\rm{\Delta }}}_{+}{x}_{i}{{\rm{\Delta }}}_{+}{x}_{j}-{{\rm{\Delta }}}_{-}{x}_{i}{{\rm{\Delta }}}_{-}{x}_{j}\right\rangle ={D}_{ij}-{D}_{ij}=0\end{array}\end{eqnarray} \tag{ 63 }$$

Also, the expression $\tfrac{\left\langle g\left(\vec{x}^{\prime\prime} \right)-g\left(\vec{x}^{\prime} \right)\right\rangle }{2{\rm{\Delta }}t}$ is the systematic derivative of $g\left(\vec{x},t\right),$ thus, we have:

$$\begin{eqnarray}&&\begin{array}{c}{D}_{c}g\left(\vec{x},t\right)=\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}+{v}_{i}^{c}\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}\end{array}\end{eqnarray} \tag{ 64 }$$

From (64) is easy to identify the systematic derivative operator as:

$$\begin{eqnarray}&&\begin{array}{c}{D}_{c}\equiv \displaystyle \frac{\partial }{\partial t}+{v}_{i}^{c}\displaystyle \frac{\partial }{\partial {x}_{i}}\end{array}\end{eqnarray} \tag{ 65 }$$

If one takes as a special case the identity function: $g\left(\vec{x},t\right)=\vec{x}\left(t\right),$ the expression is reduced to the systematic velocity in vicinity ${\mathcal B} .$

Adding (57) and (58), it results:

$$\begin{eqnarray*}&&\begin{array}{l}g\left(\vec{x}^{\prime\prime} ,t+{\rm{\Delta }}t\right)+g\left(\vec{x}^{\prime} ,t-{\rm{\Delta }}t\right)=2g\left(\vec{x},t\right)+\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}\left({x}_{i}^{{\prime\prime} }-{x}_{i}-{x}_{i}+{x}_{i}^{{\prime} }\right)\\ \,+\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\left[\left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\left({x}_{j}^{{\prime\prime} }-{x}_{j}\right)+\left({x}_{i}-{x}_{i}^{{\prime} }\right)\left({x}_{j}-{x}_{j}^{{\prime} }\right)\right]\end{array}\end{eqnarray*}$$

Rearranging, multiplying by $\tfrac{1}{2{\rm{\Delta }}t}$ and averaging, one obtains:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\left\langle g\left(\vec{x}^{\prime\prime} ,t+{\rm{\Delta }}t\right)+g\left(\vec{x}^{\prime} ,t-{\rm{\Delta }}t\right)-2g(\vec{x},t)\right\rangle }{2{\rm{\Delta }}t}=\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}\displaystyle \frac{\left\langle {{\rm{\Delta }}}_{+}{x}_{i}+{{\rm{\Delta }}}_{-}{x}_{i}\right\rangle }{2{\rm{\Delta }}t}\\ \,+\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\displaystyle \frac{\left\langle \left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\left({x}_{j}^{{\prime\prime} }-{x}_{j}\right)+\left({x}_{i}-{x}_{i}^{{\prime} }\right)\left({x}_{j}-{x}_{j}^{{\prime} }\right)\right\rangle }{2{\rm{\Delta }}t}\end{array}\end{eqnarray} \tag{ 66 }$$

The term that accompanies the first derivative with respect to the positions is used to define the diffusion velocity in vicinity ${\mathcal B} :$

$$\begin{eqnarray}&&\begin{array}{c}{u}_{i}\equiv \displaystyle \frac{\left\langle {{\rm{\Delta }}}_{+}{x}_{i}+{{\rm{\Delta }}}_{-}{x}_{i}\right\rangle }{2{\rm{\Delta }}t}\end{array}\end{eqnarray} \tag{ 67 }$$

Which is useful to measure the stirring that the state point undergoes in each subset ${\mathcal B}$ in the fluctuations space. We can back up this definition if we add the systematic velocity and the diffusion presented above, resulting in:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\left\langle \vec{x}^{\prime\prime} -\vec{x}^{\prime} \right\rangle }{2{\rm{\Delta }}t}+\displaystyle \frac{\left\langle \vec{x}^{\prime\prime} -2\vec{x}+\vec{x}^{\prime} \right\rangle }{2{\rm{\Delta }}t}=\displaystyle \frac{\left\langle \vec{x}^{\prime\prime} -\vec{x}\right\rangle }{{\rm{\Delta }}t}\end{array}\end{eqnarray} \tag{ 68 }$$

Which can be interpreted as a total velocity measuring the mean path of state points from $\vec{x}$ to $\vec{x}^{\prime\prime}$ in the time interval ${\rm{\Delta }}t.$ In that case

$$\begin{eqnarray}&&\begin{array}{c}\vec{v}={\vec{v}}_{c}+\vec{u}\end{array}\end{eqnarray} \tag{ 69 }$$

such that the cigarette smoke phenomenon can be described in terms of two velocities, one of translation and other of swirling for each of the ${\mathcal B}$ vicinities in the fluctuation space.

With these conceptual tools, one can identify the operators that allow the calculation of the velocities mentioned in this work. Combining results, we find the next expression:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\left\langle g\left(\vec{x}^{\prime\prime} ,t+{\rm{\Delta }}t\right)+g\left(\vec{x}^{\prime} ,t-{\rm{\Delta }}t\right)-2g\left(\vec{x},t\right)\right\rangle }{2{\rm{\Delta }}t}+{u}_{i}\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\end{array}\end{eqnarray} \tag{ 70 }$$

where

$$\begin{eqnarray}&&\begin{array}{c}2{D}_{ij}=\displaystyle \frac{\left\langle {{\rm{\Delta }}}_{-}{x}_{i}{{\rm{\Delta }}}_{-}{x}_{j}+{{\rm{\Delta }}}_{-}{x}_{i}{{\rm{\Delta }}}_{-}{x}_{j}\right\rangle }{2{\rm{\Delta }}t}\end{array}\end{eqnarray} \tag{ 71 }$$

With the elements considered up until now, we identify the left-hand side as the stochastic derivative, or of diffusion, of the function $g(\vec{x},t)$ and identify the operator of the stochastic derivative as

$$\begin{eqnarray}&&\begin{array}{c}{D}_{s}={u}_{i}\displaystyle \frac{\partial }{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\end{array}\end{eqnarray} \tag{ 72 }$$

such that

$$\begin{eqnarray}&&\begin{array}{c}{D}_{s}g\left(\vec{x},t\right)={u}_{i}\end{array}\end{eqnarray} \tag{ 73 }$$

To study the exit velocity, ${\vec{v}}^{e},$ we study the displacement from $\vec{x}$ to $\vec{x}^{\prime\prime}$ in a ${\rm{\Delta }}t$ time interval. For this purpose, we reuse the Taylor expansion given in (57) and rearrange it as:

$$\begin{eqnarray*}&&\begin{array}{l}g\left(\vec{x}^{\prime\prime} ,t+{\rm{\Delta }}t\right)-g\left(\vec{x},t\right)=\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}{\rm{\Delta }}t+\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}\left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\\ \,+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\left({x}_{j}^{{\prime\prime} }-{x}_{j}\right)+\ldots \end{array}\end{eqnarray*}$$

Multiplying by $\tfrac{1}{{\rm{\Delta }}t}$ and calculating the mean

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}\displaystyle \frac{\left\langle g\left(\vec{x}^{\prime\prime} ,t+{\rm{\Delta }}t\right)-g\left(\vec{x},t\right)\right\rangle }{{\rm{\Delta }}t}=\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}+\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}\displaystyle \frac{\left\langle {x}_{i}^{{\prime\prime} }-{x}_{i}\right\rangle }{{\rm{\Delta }}t}\\ \,+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\displaystyle \frac{\left\langle \left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\left({x}_{j}^{{\prime\prime} }-{x}_{j}\right)\right\rangle }{{\rm{\Delta }}t}+\ldots \end{array}\end{array}\end{eqnarray} \tag{ 74 }$$

Defining the $i$-th component of the exit velocity as follows

$$\begin{eqnarray}&&\begin{array}{c}{v}_{i}^{e}\equiv \displaystyle \frac{\left\langle {x}_{i}^{{\prime\prime} }-{x}_{i}\right\rangle }{{\rm{\Delta }}t}\end{array}\end{eqnarray} \tag{ 75 }$$

and identifying that

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{1}{2}\displaystyle \frac{\left\langle \left({x}_{i}^{{\prime\prime} }-{x}_{i}\right)\left({x}_{j}^{{\prime\prime} }-{x}_{j}\right)\right\rangle }{{\rm{\Delta }}t}=\displaystyle \frac{1}{2{\rm{\Delta }}t}\left\langle {{\rm{\Delta }}}_{+}{x}_{i}{{\rm{\Delta }}}_{+}{x}_{j}\right\rangle ={D}_{ij}\end{array}\end{eqnarray} \tag{ 76 }$$

thus, we have

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\left\langle g\left(\vec{x}^{\prime\prime} ,t+{\rm{\Delta }}t\right)-g\left(\vec{x},t\right)\right\rangle }{{\rm{\Delta }}t}=\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}+{v}_{i}^{e}\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\end{array}\end{eqnarray} \tag{ 77 }$$

Defining the left-hand side as a forward derivative of function $g\left(\vec{x},t\right)$ and introduce the forward operator as:

$$\begin{eqnarray}&&\begin{array}{c}{D}^{e}\equiv \displaystyle \frac{\partial }{\partial t}+{v}_{i}^{e}\displaystyle \frac{\partial }{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\end{array}\end{eqnarray} \tag{ 78 }$$

Such that we represent the forward derivative as ${D}^{e}g\left(\vec{x},t\right).$ In the case when $g\left(\vec{x},t\right)={x}_{i}$ one obtains the $i$-th component of the exit velocity in vicinity ${\mathcal B} .$

$$\begin{eqnarray}&&\begin{array}{c}{D}^{e}{x}_{i}={v}_{i}^{e}\end{array}\end{eqnarray} \tag{ 79 }$$

Finally, one can relate the entry velocity with the translation of the state point from $\vec{x}^{\prime}$ to $\vec{x}$ in time interval ${\rm{\Delta }}t$ (see figure 5). Now consider the Taylor series expansion:

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}g\left(\vec{x}^{\prime} ,t-{\rm{\Delta }}t\right)=g\left(\vec{x},t\right)-\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}{\rm{\Delta }}t-\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}\left({x}_{i}-{x}_{i}^{{\prime} }\right)\\ \,+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\left({x}_{i}-{x}_{i}^{{\prime} }\right)\left({x}_{j}-{x}_{j}^{{\prime} }\right)+\ldots \end{array}\end{array}\end{eqnarray} \tag{ 80 }$$

Rearranging, multiplying by $\tfrac{1}{{\rm{\Delta }}t}$ and calculating the mean:

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}\displaystyle \frac{\left\langle g\left(\vec{x}^{\prime} ,t-{\rm{\Delta }}t\right)-g\left(\vec{x},t\right)\right\rangle }{{\rm{\Delta }}t}=-\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}-\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}\displaystyle \frac{\left\langle {x}_{i}-{x}_{i}^{{\prime} }\right\rangle }{{\rm{\Delta }}t}\\ \,+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\displaystyle \frac{\left\langle \left({x}_{i}-{x}_{i}^{{\prime} }\right)\left({x}_{j}-{x}_{j}^{{\prime} }\right)\right\rangle }{{\rm{\Delta }}t}\end{array}\end{array}\end{eqnarray} \tag{ 81 }$$

Defining the ith component of the entry velocity in vicinity ${\mathcal B}$ as:

$$\begin{eqnarray}&&\begin{array}{c}{v}_{i}^{a}\equiv \displaystyle \frac{\left\langle {x}_{i}-{x}_{i}^{{\prime} }\right\rangle }{{\rm{\Delta }}t}\end{array}\end{eqnarray} \tag{ 82 }$$

and identifying

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{1}{2}\displaystyle \frac{\left\langle \left({x}_{i}-{x}_{i}^{{\prime} }\right)\left({x}_{j}-{x}_{j}^{{\prime} }\right)\right\rangle }{{\rm{\Delta }}t}={D}_{ij}\end{array}\end{eqnarray} \tag{ 83 }$$

It is possible to rewrite (81) as:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\left\langle g\left({\vec{x}}^{^{\prime} },t-{\rm{\Delta }}t\right)-g\left(\vec{x},t\right)\right\rangle }{{\rm{\Delta }}t}=-\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}-{v}_{i}^{a}\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\end{array}\end{eqnarray} \tag{ 84 }$$

Defining the backwards derivative operator as

$$\begin{eqnarray}&&\begin{array}{c}{D}_{a}\equiv -\left(-\displaystyle \frac{\partial }{\partial t}-{v}_{i}^{a}\displaystyle \frac{\partial }{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\right)\end{array}\end{eqnarray} \tag{ 85 }$$

And rewriting (84) as

$$\begin{eqnarray}&&\begin{array}{c}{D}_{a}g\left(\vec{x},t\right)=-\left(-\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial t}-{v}_{i}^{a}\displaystyle \frac{\partial g\left(\vec{x},t\right)}{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}g\left(\vec{x},t\right)}{\partial {x}_{i}\partial {x}_{j}}\right)\end{array}\end{eqnarray} \tag{ 86 }$$

Once again, if $g\left(\vec{x},t\right)={x}_{i}$ one can obtain

$$\begin{eqnarray}&&\begin{array}{c}{D}_{a}{x}_{i}={v}_{i}^{a}\end{array}\end{eqnarray} \tag{ 87 }$$

which is the ith component of the entry velocity in vicinity ${\mathcal B} .$

Combining operators

Passing the operators through an algebraic process, one can obtain the following expressions:

$$\begin{eqnarray}&&\begin{array}{c}{D}^{e}\equiv \displaystyle \frac{\partial }{\partial t}+{v}_{i}^{e}\displaystyle \frac{\partial }{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\end{array}\end{eqnarray} \tag{ 88 }$$

$$\begin{eqnarray}&&\begin{array}{c}{D}_{a}\equiv \displaystyle \frac{\partial }{\partial t}+{v}_{i}^{a}\displaystyle \frac{\partial }{\partial {x}_{i}}-{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\end{array}\end{eqnarray} \tag{ 89 }$$

If we were to add (88) and (89), then multiply by $\tfrac{1}{2}$

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{1}{2}\left({D}^{e}+{D}_{a}\right)=\displaystyle \frac{\partial }{\partial t}+\displaystyle \frac{{v}_{i}^{e}+{v}_{i}^{a}}{2}\displaystyle \frac{\partial }{\partial {x}_{i}}\end{array}\end{eqnarray} \tag{ 90 }$$

Defining

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{{v}_{i}^{e}+{v}_{i}^{a}}{2}\equiv {v}_{i}^{c}\end{array}\end{eqnarray} \tag{ 91 }$$

we have

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{1}{2}\left({D}^{e}+{D}_{a}\right)=\displaystyle \frac{\partial }{\partial t}+\displaystyle \frac{{v}_{i}^{e}+{v}_{i}^{a}}{2}\displaystyle \frac{\partial }{\partial {x}_{i}}=\displaystyle \frac{\partial }{\partial t}+{v}_{i}^{c}\displaystyle \frac{\partial }{\partial {x}_{i}}={D}^{c}\end{array}\end{eqnarray} \tag{ 92 }$$

Now, calculating the difference between (88) and (89) and multiplying by $\tfrac{1}{2}$

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{1}{2}\left({D}^{e}-{D}_{a}\right)=\displaystyle \frac{{v}_{i}^{e}-{v}_{i}^{a}}{2}\displaystyle \frac{\partial }{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\end{array}\end{eqnarray} \tag{ 93 }$$

such that it is convenient to define:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{{v}_{i}^{e}-{v}_{i}^{a}}{2}\equiv {u}_{i}\end{array}\end{eqnarray} \tag{ 94 }$$

So (93) can be rewritten as:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{1}{2}\left({D}^{e}-{D}_{a}\right)={u}_{i}\displaystyle \frac{\partial }{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}={D}^{s}\end{array}\end{eqnarray} \tag{ 95 }$$

The stochastic velocities provide a description at the level of vicinities such that it is possible to calculate a variety of means previously mentioned. Table 2 summarizes the previous results:

Table 2. Stochastic velocities. Notation and associated operator.

Velocity	Notation	Operator
Access entry	${\vec{v}}^{a}$	${D}^{e}\vec{x}\left(t\right)=\left(\displaystyle \frac{\partial }{\partial t}+{v}_{i}^{e}\displaystyle \frac{\partial }{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\right)\vec{x}\left(t\right)$
Exit	${\vec{v}}^{e}$	${D}^{a}\vec{x}\left(t\right)=\left(\displaystyle \frac{\partial }{\partial t}+{v}_{i}^{a}\displaystyle \frac{\partial }{\partial {x}_{i}}-{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\right)\vec{x}\left(t\right)$
Systematic	${\vec{v}}^{c}$	${D}^{c}\vec{x}\left(t\right)=\left(\displaystyle \frac{\partial }{\partial t}+{v}_{i}^{c}\displaystyle \frac{\partial }{\partial {x}_{i}}\right)\vec{x}\left(t\right)$
Diffusion	$\vec{u}$	${D}^{s}\vec{x}\left(t\right)=\left({u}_{i}\displaystyle \frac{\partial }{\partial {x}_{i}}+{D}_{ij}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}\right)\vec{x}\left(t\right)$

In their current form, their usefulness is not very clear. In the following section we will see a version of these that allows us to study the quasi-local conduct of the state points in the fluctuation space. Figures 5 and 7 lead to a description where the idea of instantaneous velocities, used regularly in the context of classical mechanics, has to be abandoned. In the topic under study there is a necessity to associate the concepts of velocities to vicinities that are small enough, but without reducing their size to an infinitesimal area, as is the standard in differential calculus.

The stochastic velocities have a practical application in time dependent Ornstein-Uhlenbeck stochastic processes because they can be written in terms of the convection, diffusion and self-correlation matrices of these processes. To make the understanding of the conduct of the state point in the fluctuation space more accessible, in this section we obtain the form of these velocities for the system under consideration.

The time dependent Ornstein-Uhlenbeck processes are normally distributed, $P\left(\vec{q},t\right),$ in which their means and self-correlation function change with time. The forward Fokker-Planck equation (FPE) that satisfies $P\left(\vec{q},t\right)$ is written as follows:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\partial P\left(\vec{q},t\right)}{dt}=\displaystyle \frac{\partial \left[{F}_{\mu }\left(\vec{q},t\right)P\left(\vec{q},t\right)\right]}{\partial {q}_{\mu }}+\displaystyle \frac{{\partial }^{2}\left[{D}_{\mu \nu }\left(t\right)P\left(\vec{q},t\right)\right]}{\partial {q}_{\mu }\partial {q}_{\nu }}\end{array}\end{eqnarray} \tag{ 96 }$$

The factor of $\tfrac{1}{2}$ frequently used when writing a FPE has been absorbed by ${D}_{\mu \nu },$ and the repeated indexes run from $1$ to $p,$ with $p$ denoting the number of degrees of liberty of the system. The flux term, ${F}_{\mu },$ is linear in the noise $\vec{q},$ as shown below:

$$\begin{eqnarray}&&\begin{array}{c}{F}_{\mu }\left(\vec{q},t\right)=-{L}_{\mu \nu }{q}_{\nu }\end{array}\end{eqnarray} \tag{ 97 }$$

where ${L}_{\mu \nu }$ is called the convection matrix.

The probability distribution that satisfies the forward FPE is given in (98)

$$\begin{eqnarray}&&\begin{array}{c}P\left(\vec{q},t\right)=\displaystyle \frac{1}{{\left(2\pi \right)}^{p/2}\sqrt{Det\left({\rm{\Xi }}\left(t\right)\right)}}{e}^{-\displaystyle \frac{1}{2}\left({q}_{\mu }-\left\langle {q}_{\mu }\left(t\right)\right\rangle \right){\left[{{\rm{\Xi }}}^{-1}\left(t\right)\right]}_{\mu \nu }\left({q}_{\nu }-\left\langle {q}_{\nu }\left(t\right)\right\rangle \right)}\end{array}\end{eqnarray} \tag{ 98 }$$

We can identify the exit velocity ${\vec{v}}_{e}$ with the flux term, thus we have:

$$\begin{eqnarray}&&\begin{array}{c}{v}_{\mu }^{e}=-{L}_{\mu \nu }{q}_{\nu }\end{array}\end{eqnarray} \tag{ 99 }$$

so, we can rewrite the forward FPE as

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\partial P\left(\vec{q},t\right)}{dt}=-{v}_{\mu }^{e}\displaystyle \frac{\partial P\left(\vec{q},t\right)}{\partial {q}_{\mu }}+{D}_{\mu \nu }\left(t\right)\displaystyle \frac{{\partial }^{2}P\left(\vec{q},t\right)}{\partial {q}_{\mu }\partial {q}_{\nu }}\end{array}\end{eqnarray} \tag{ 100 }$$

On the other hand, the backward FPE that corresponds to this process takes the following form:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\partial P\left(\vec{q},t\right)}{dt}=-{F}_{\mu }\left(\vec{q},t\right)\displaystyle \frac{\partial P\left(\vec{q},t\right)}{\partial {q}_{\mu }}-{D}_{\mu \nu }\left(t\right)\displaystyle \frac{{\partial }^{2}P\left(\vec{q},t\right)}{\partial {q}_{\mu }\partial {q}_{\nu }}\end{array}\end{eqnarray} \tag{ 101 }$$

In a similar fashion, the backward FPE is related to the operator ${D}_{a}$ and with the entry velocity, resulting in

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\partial P\left(\vec{q},t\right)}{dt}=-{v}_{\mu }^{a}\displaystyle \frac{\partial P\left(\vec{q},t\right)}{\partial {q}_{\mu }}-{D}_{\mu \nu }\left(t\right)\displaystyle \frac{{\partial }^{2}P\left(\vec{q},t\right)}{\partial {q}_{\mu }\partial {q}_{\nu }}\end{array}\end{eqnarray} \tag{ 102 }$$

There is an analytical solution when the diffusion velocity has zero divergence. Calculating the difference between (100) and (102) and multiplying by $\tfrac{1}{2}$

$$\begin{eqnarray}&&\begin{array}{c}0=\displaystyle \frac{\partial }{\partial {q}_{\mu }}\left[{u}_{\mu }P\left(\vec{q},t\right)-{D}_{\mu \nu }\left(t\right)\displaystyle \frac{\partial P\left(\vec{q},t\right)}{\partial {q}_{\nu }}\right]\end{array}\end{eqnarray} \tag{ 103 }$$

The term between brackets is a magnitude with divergence of zero:

$$\begin{eqnarray}&&\begin{array}{c}{J}_{\mu }={u}_{\mu }P\left(\vec{q},t\right)-{D}_{\mu \nu }\left(t\right)\displaystyle \frac{\partial P\left(\vec{q},t\right)}{\partial {q}_{\nu }}\end{array}\end{eqnarray} \tag{ 104 }$$

thus, it can be interpreted as a magnitude that is conserved:

$$\begin{eqnarray}&&\begin{array}{c}{J}_{\mu }={u}_{\mu }P\left(\vec{q},t\right)-{D}_{\mu \nu }\left(t\right)\displaystyle \frac{\partial P\left(\vec{q},t\right)}{\partial {q}_{\nu }}=C\end{array}\end{eqnarray} \tag{ 105 }$$

where $C$ is a constant that can be taken as equal to zero, then ${u}_{\mu }$ is given in terms of $P\left(\vec{q},t\right)$ and ${D}_{\mu \nu }:$

$$\begin{eqnarray}&&\begin{array}{c}{u}_{\mu }={D}_{\mu \nu }\left(t\right)\displaystyle \frac{1}{P\left(\vec{q},t\right)}\displaystyle \frac{\partial P\left(\vec{q},t\right)}{\partial {q}_{\nu }}={D}_{\mu \nu }\left(t\right)\displaystyle \frac{\partial \,\mathrm{ln}\left|P\left(\vec{q},t\right)\right|}{\partial {q}_{\nu }}\end{array}\end{eqnarray} \tag{ 106 }$$

From above it results that the diffusion velocity can be obtained if $P\left(\vec{q},t\right)$ is known. This is the case for time dependent Ornstein-Uhlenbeck processes.

Substituting (98) in (106) and using the short notation of: $G\left(t\right)=\tfrac{1}{2}{y}_{\alpha }\left(t\right){\left[{{\rm{\Xi }}}^{-1}\left(t\right)\right]}_{\alpha \beta }{y}_{\beta }\left(t\right),$ one obtains

$$\begin{eqnarray}&&\begin{array}{c}{u}_{\mu }={D}_{\mu \nu }\left(t\right)\displaystyle \frac{\partial }{\partial {q}_{\nu }}\,\mathrm{ln}\left\{\displaystyle \frac{1}{{\left(2\pi \right)}^{p/2}\sqrt{Det\left({\rm{\Xi }}\left(t\right)\right)}}{e}^{-G\left(t\right)}\right\}\end{array}\end{eqnarray} \tag{ 107 }$$

where we have written $y\left(t\right)={q}_{\mu }\left(t\right)-\left\langle q\_\mu \left(t\right)\right\rangle$ to simplify notation. Working (107) some more results in:

$$\begin{eqnarray}&&\begin{array}{c}{u}_{\mu }=-{D}_{\mu \nu }\left(t\right)\displaystyle \frac{\partial G\left(t\right)}{\partial {q}_{\nu }}\end{array}\end{eqnarray} \tag{ 108 }$$

Calculating the derivative and using that ${\left({{\rm{\Xi }}}^{-1}\right)}_{\alpha \beta }$ is symmetric results that

$$\begin{eqnarray}&&\begin{array}{c}{u}_{\mu }=-{D}_{\mu \nu }\left(t\right){\left[{{\rm{\Xi }}}^{-1}\left(t\right)\right]}_{\nu \beta }\left[{q}_{\beta }\left(t\right)-\left\langle {q}_{\beta }\left(t\right)\right\rangle \right]\end{array}\end{eqnarray} \tag{ 109 }$$

From previous results we have the next relations for the velocities:

$$\begin{eqnarray}&&\begin{array}{c}{\vec{v}}^{c}=\displaystyle \frac{1}{2}\left({\vec{v}}^{e}+{\vec{v}}^{a}\right)\end{array}\end{eqnarray} \tag{ 110 }$$

$$\begin{eqnarray}&&\begin{array}{c}\vec{u}=\displaystyle \frac{1}{2}\left({\vec{v}}^{e}-{\vec{v}}^{a}\right)\end{array}\end{eqnarray} \tag{ 111 }$$

Relation (111) leads to

$$\begin{eqnarray*}&&{\vec{v}}^{a}=-2\vec{u}+{\vec{v}}^{e}\end{eqnarray*}$$

Substituting this result in (110) gives:

$$\begin{eqnarray*}&&{\vec{v}}^{c}={\vec{v}}^{e}-\vec{u}\end{eqnarray*}$$

so, the components of the entry velocity have the form:

$$\begin{eqnarray}&&\begin{array}{c}{v}_{\mu }^{a}=-2{u}_{\mu }+{v}_{\mu }^{e}=2{D}_{\mu \nu }{\left[{{\rm{\Xi }}}^{-1}\left(t\right)\right]}_{\nu \beta }\left[{q}_{\beta }\left(t\right)-\left\langle {q}_{\beta }\left(t\right)\right\rangle \right]-{L}_{\mu \nu }{q}_{\nu }\end{array}\end{eqnarray} \tag{ 112 }$$

and the components of the systematic velocity are:

$$\begin{eqnarray}&&\begin{array}{c}{v}_{\mu }^{c}=-{L}_{\mu \nu }{q}_{\nu }+{D}_{\mu \nu }\left(t\right){\left[{{\rm{\Xi }}}^{-1}\left(t\right)\right]}_{\nu \beta }\left[{q}_{\beta }\left(t\right)-\left\langle {q}_{\beta }\left(t\right)\right\rangle \right]\end{array}\end{eqnarray} \tag{ 113 }$$

With this, the set of stochastic velocities in terms of the probability density is now complete. The results above can be written in matrix notation. Table 3 displays the four stochastic velocities for time dependent Ornstein-Uhlenbeck processes:

Table 3. Stochastic velocity operators for Ornstein-Uhlenbeck processes.

Velocity	Notation	Operator
entry	${\vec{v}}^{a}$	$2D{{\rm{\Xi }}}^{-1}\left(t\right)\left(\vec{q}-\left\langle \vec{q}\right\rangle \right)-L\vec{q}$
Exit	${\vec{v}}^{e}$	${\vec{v}}^{e}=-L\vec{q}$
Systematic	${\vec{v}}^{c}$	$D{{\rm{\Xi }}}^{-1}\left(t\right)\left(\vec{q}-\left\langle \vec{q}\right\rangle \right)-L\vec{q}$
Diffusion	$\vec{u}$	$-D{{\rm{\Xi }}}^{-1}\left(t\right)\left(\vec{q}-\left\langle \vec{q}\right\rangle \right)$

Table 4. Values of the mass and number of molecules for enzyme and substrate.

${M}_{enz}=4.981616763\times {10}^{-23}\,{\rm{kg}}$	${M}_{sus}=5.73314\times {10}^{-25}\,{\rm{kg}}$
${N}_{enz}=100$	${N}_{sus}=4999$

The diffusion velocity and the systematic velocity will be used below to describe the behavior of a state point in each vicinity ${\mathcal B}$ in the fluctuation space.

The analysis of the simulated catalysis reaction allows the determination of when thermodynamic equilibrium has been reached. From expression (37), it is evident that the determinant can be used for that purpose. The state of equilibrium is reached if $\det \left({\rm{\Xi }}\left({t}_{2}\right)\right)=\det \left({\rm{\Xi }}\left({t}_{1}\right)\right),$ where ${t}_{1}\lt {t}_{2}.$ In numerical calculation it is possible to define a tolerance

$$\begin{eqnarray*}&&{\rm{tol}}=\displaystyle \frac{{\rm{\det }}({\rm{\Xi }}\left({t}_{2}\right))}{{\rm{\det }}({\rm{\Xi }}\left({t}_{1}\right))}\lt {10}^{-5}\end{eqnarray*}$$

This is the regime called equilibrium state, and it corresponds to the moment when substrate and enzyme-substrate complex have been depleted. It is achieved at $t\geqslant 0.013.$

The quasi-stationary process is studied through the numerical solution of equations (32). The expressions for the time evolution of the substrate and enzyme-substrate complex presented in table 6 in section Entropy of Fluctuations are utilized to calculate the autocorrelation functions shown in figure 8 below

Zoom In Zoom Out Reset image size

Figure 8(a) exhibits the autocorrelation of the fluctuation of the substrate, figure 8(b) is the graph of the autocorrelation of the fluctuation of the enzyme-substrate complex, and figure 8(c) the correlation of both. The numerical analysis could be performed up to $t=0.015,$ but the description would no longer correspond to the state under study, the quasi-stationary state; instead, it would be a state with only the remnants of the noise of the substrate and complex concentrations, which is of little practical interest in this work.

The probability density can be calculated with expression (27). Figure 9 shows its initial form at $t=0.0001,$ and its final form at $t=0.0158$ once it reaches the end of the quasi-stationary state. At first glance the differences between the gaussian distributions at its initial and final states are not apparent. But a more careful observation reveals that there has been a clockwise rotation of approximately $-4^\circ ,$ as shown in figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An analytical way of detecting a change is by calculating the difference between probability densities at ${t}_{f}=0.0158$ and ${t}_{i}=0.0001,$ $P\left(\vec{q},{t}_{f}\right)-P\left(\vec{q},{t}_{i}\right).$ The result can be seen in figure 11, where there is a region at the center with higher probability, while above and below there is a region with lower probability. What this means is that, during the quasi-stationary state of the catalytic process, the probability density shifts towards the center.

Zoom In Zoom Out Reset image size

This is an expected outcome once proper attention is paid to the graphs of ${{\rm{\Xi }}}_{11}\left(t\right)$ and ${{\rm{\Xi }}}_{22}\left(t\right),$ where the former increases as the latter decreases with time. Further details on information contained within the quasi-station state, that is relevant to the understanding of the biochemical process, will be addressed in a later section titled Entropy of Fluctuations.

Before going to the next section, where we address the topic of the entropy in this model, let us focus on the total stochastic velocity, which results from the sum of the systematic velocity and the diffusion velocity. Figure 12 shows a comparison of the initial state, at $t=0.0001,$ and a state close to equilibrium, at $t=0.014,$ of the total stochastic velocity.

Zoom In Zoom Out Reset image size

The total stochastic velocity, shown in figure 12, is a two-dimensional field that plots the tendency of the state points $\vec{q}$ to move towards a region in the center, which intuitively coincides with the plot of the probability density. At the start of the reaction ($t=0.0001$) it is mainly the points from the second and fourth quadrant that move at a greater velocity, contributing the most in keeping the saturation on the center. In contrast, when very close to the equilibrium ($t=0.014$) the roles have reversed, and it is now the first and third quadrants responsible of keeping the saturation in the center.

We calculate the curl of the total stochastic velocity to further explore this behavior, as shown in figure 13. From a geometrical viewpoint, the tendency of the field to rotate suggests that the averaged transitions performed by the state points, $\vec{q},$ within each region ${\mathcal B}$ occur due to the lack of balance between the displacements $\pm {q}_{x},$ which is something that is also reflected by the averaged $\pm {q}_{y}.$ This conduct causes the semi-axes of symmetry of the probability distribution to not remain static. The magnitude of this phenomenon diminishes with time, corresponding to the end of the quasi-stationary stage.

Zoom In Zoom Out Reset image size

The second law of thermodynamics establishes that, in the absence of external work being done on a system, entropy follows the inequality shown in (116):

$$\begin{eqnarray}&&\begin{array}{c}{\rm{\Delta }}S\geqslant 0\end{array}\end{eqnarray} \tag{ 116 }$$

Such that the equal sign is present once the state of equilibrium is reached.

Instead, entropy can diminish if energy is applied to a system through appropriate processes. We see that this is the case in enzymatic catalysis that are studied with the Michaelis-Menten model.

From the definition of Markov processes [34], it is clear that, in mathematical terms, the decrease of entropy through time is not forbidden. For a physical system with microscopic states numbered with $n,$ and its probability is denoted by ${P}_{n}\left(t\right),$ its entropy can be expressed as seen in (117):

$$\begin{eqnarray}&&\begin{array}{c}S\left(t\right)=-\displaystyle \displaystyle \sum _{n}{P}_{n}\left(t\right)\mathrm{ln}\,{{\rm{P}}}_{{\rm{n}}}\left({\rm{t}}\right)\end{array}\end{eqnarray} \tag{ 117 }$$

Its rate of change is given by (118):

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{dS\left(t\right)}{dt}=-\displaystyle \displaystyle \sum _{n}\displaystyle \frac{\partial {P}_{n}\left(t\right)}{\partial t}\left[\mathrm{ln}\,{P}_{n}\left(t\right)+1\right]\end{array}\end{eqnarray} \tag{ 118 }$$

From $0\leqslant {P}_{n}\left(t\right)\leqslant 1,$ results that $-\infty \lt \,\mathrm{ln}\,{P}_{n}\left(t\right)\leqslant 0.$ Therefore, the inequality $\tfrac{dS\left(t\right)}{dt}\lt 0$ can be fulfilled if either of the conditions (119) or (120) are followed.

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\partial {P}_{n}\left(t\right)}{\partial t}\lt 0\,{\rm{and}}\,\mathrm{ln}\,{P}_{n}\left(t\right)\lt -1\end{array}\end{eqnarray} \tag{ 119 }$$

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\partial {P}_{n}\left(t\right)}{\partial t}\gt 0\,{\rm{and}}\,\mathrm{ln}\,{P}_{n}\left(t\right)\gt -1\end{array}\end{eqnarray} \tag{ 120 }$$

In this section we study the entropy in the time interval $t\leqslant 0$ up to equation (135). From that point on, the contribution of the entropy of fluctuations is added to the Michaelis-Menten model of the penicillin hydrolysis by the beta-lactamase. This is the case when the initial state of equilibrium is broken, then begins a process where the most important stage is the quasi-stationary state of the reaction, that reaches a state of equilibrium once the substrate has been exhausted.

The entropy of a system in statistical physics is calculated by $S={k}_{B}{\left(\tfrac{\partial T\mathrm{ln}Q}{\partial T}\right)}_{N,V},$ where $Q$ is the partition function. But as stated before, we will suppose that the system behaves like an ideal gas system. Therefore, for an ideal gas of $N$ particles of mass $m,$ the expression for the entropy is given by:

$$\begin{eqnarray*}&&{S}_{ideal}^{\left(i\right)}={N}_{i}{k}_{B}\,\mathrm{ln}\left[{\left(\displaystyle \frac{2\pi {m}_{i}{k}_{B}T}{{h}^{2}}\right)}^{3/2}\displaystyle \frac{{q}^{5/2}}{{c}_{i}}\right]\end{eqnarray*}$$

Where ${c}_{i}={c}_{1},\,{c}_{2},$ with ${c}_{1}=\tfrac{{N}_{S}}{{\rm{\Omega }}},\,{c}_{2}=\tfrac{{N}_{E}}{{\rm{\Omega }}},$ and ${N}_{S}$ the number of substrate particles, ${N}_{E}$ the number of enzyme particles, and ${\rm{\Omega }}$ is the initial number of substrate and enzymes in the system. In this expression ${m}_{i}$ denotes the masses, ${m}_{1}$ is the mass of a substrate particle and ${m}_{2}$ is the mass of the enzyme particle.

Translation entropy

For the substrate, let us consider as an example the antibiotic rifampicin, which is part of the family of penicillin. The enzyme considered is the beta-lactamase. The data used to model these molecules is shown in table 4:

With a temperature of $T=36.5\,^\circ C=309.65\,{\rm{K}},$ the translation entropy of substrate and enzyme are given in (121) and (122), respectively.

$$\begin{eqnarray}&&\begin{array}{c}{S}_{tr}^{\left(sus\right)}=80.3875\,{N}_{sus}{k}_{B}\end{array}\end{eqnarray} \tag{ 121 }$$

$$\begin{eqnarray}&&\begin{array}{c}{S}_{tr}^{\left(enz\right)}=1.82029\,\eta {N}_{sus}{k}_{B}\end{array}\end{eqnarray} \tag{ 122 }$$

Where $\eta =\tfrac{{N}_{enz}}{{N}_{sus}}=\tfrac{100}{4999}=0.02.$ The total translation entropy is given by (123):

$$\begin{eqnarray}&&\begin{array}{c}{S}_{tr}={S}_{tr}^{\left(sus\right)}+{S}_{tr}^{\left(enz\right)}={N}_{sus}{k}_{B}\left[80.3875+\eta 1.82029\right]\end{array}\end{eqnarray} \tag{ 123 }$$

Mixing entropy

This contribution to the total entropy exists due to the existence of two gasses mixing inside a volume. The mixing entropy can be expressed as:

$$\begin{eqnarray*}&&{S}_{mix}=-{\rm{\Omega }}{k}_{B}\left({c}_{sus}\,\mathrm{ln}\,{c}_{sus}+{c}_{enz}\,\mathrm{ln}\,{c}_{enz}\right)\end{eqnarray*}$$

Where ${c}_{sus}=\tfrac{{N}_{sus}}{{\rm{\Omega }}}=\tfrac{4999}{4999+100}=\tfrac{4999}{5099}=0.98039$ and ${c}_{enz}=\tfrac{{N}_{enz}}{{\rm{\Omega }}}=1.9612\times {10}^{-2}.$ The estimate of the mixing entropy is given by (124):

$$\begin{eqnarray}&&\begin{array}{c}{S}_{mix}=0.0984549\left(1+\eta \right){N}_{sus}{k}_{B}\end{array}\end{eqnarray} \tag{ 124 }$$

Vibrational entropy

Here we study the contribution to the entropy by the vibrations of enzyme and substrate molecules. We begin with the substrate.

Vibrational entropy of the substrate

The rifampicin molecule is relatively small when compared to that of the beta-lactamase. For this analysis we suppose that the molecules can be considered as oscillating nuclei that can be detected by IR spectroscopy, therefore we model the substrate as $p$ independent harmonic oscillators, where ${\nu }_{k}$ denotes the frequencies, with $k=1,2,\ldots ,p.$ The contribution to the entropy due to the vibrational degrees of freedom of the substrate can be calculated with:

$$\begin{eqnarray*}&&{S}_{v}^{\left(sus\right)}={k}_{B}\left\{\displaystyle \frac{\partial }{\partial T}\left[T\,\mathrm{ln}\left(\displaystyle \prod _{p}^{k=1}{\left(\displaystyle \frac{1}{1-{e}^{-\beta h{\nu }_{k}}}\right)}^{{N}_{sus}}\right)\right]\right\}\end{eqnarray*}$$

This results in (125):

$$\begin{eqnarray}&&\begin{array}{c}{S}_{v}^{\left(sus\right)}={N}_{sus}{k}_{B}\left[-\displaystyle \sum _{k=1}^{p}\mathrm{ln}\left(1-{e}^{-\displaystyle \frac{h{\nu }_{k}}{{k}_{B}T}}\right)+\displaystyle \sum _{p}^{k=1}\displaystyle \frac{\displaystyle \frac{h{\nu }_{k}}{{k}_{B}T}}{{e}^{-\displaystyle \frac{h{\nu }_{k}}{{k}_{B}T}}-1}\right]\end{array}\end{eqnarray} \tag{ 125 }$$

Ivashchenko [35] provides a set of wave numbers that, when transformed to frequencies, provide the following values:

$$\begin{eqnarray*}&&v=\left\{3.71,\,4.12,\,4.30,\,4.36,\,4.60,\,4.68,\,4.87,\,4.95,\,5.16,\,8.77,\,9.14,\,10.04\right\}\times {10}^{13}\,{\rm{Hz}}\end{eqnarray*}$$

For a temperature of $T=36.5\,^\circ C=309.65\,\,{\rm{K}}$ the vibrational entropy contribution of the substrate is shown in (126):

$$\begin{eqnarray}&&\begin{array}{c}{S}_{v}^{\left(sus\right)}=12.0203\,{N}_{sus}{k}_{B}\end{array}\end{eqnarray} \tag{ 126 }$$

Vibrational entropy of the enzyme

The enzyme is a sizeable molecule, compared to the substrate, but its dimensions are still within the order of nanometers. For it we use the oscillatory theory of very small solids with the added correction necessary for such scales. The shape also plays an important role, but we will suppose it is a sphere of diameter between 3 to 7 nm, with a homogenously distributed mass.

Bu-Xuan Wang [36] argues that the theory of specific heats of Einstein can be applied to nanoparticles, thus, we use the partition function of solids [37] given as:

$$\begin{eqnarray*}&&{Q}_{v}^{\left(enz\right)}=-\displaystyle {\int }_{0}^{\infty }\left[\mathrm{ln}\left(1-{e}^{-\beta h\nu }\right)+\displaystyle \frac{1}{2}\beta h\nu \right]g\left(\nu \right)d\nu \end{eqnarray*}$$

Using the hypothesis by Einstein of $g\left(\nu \right)=3{N}_{enz}\delta \left(\nu -{\nu }_{E}\right),$ with ${\nu }_{E}$ being the Einstein frequency, and using the empirical adjustment parameter known as the Einstein temperature, ${{\rm{\Theta }}}_{E}=\tfrac{h{\nu }_{E}}{{k}_{B}},$ the expression for the vibrational entropy of the enzyme can be obtained, as shown in (127):

$$\begin{eqnarray}&&\begin{array}{l}{S}_{v}^{\left(enz\right)}=-3{k}_{B}{N}_{enz}\left\{\displaystyle \frac{\partial }{\partial T}\left[T\,\mathrm{ln}\left(1-{e}^{-{{\rm{\Theta }}}_{E}/T}\right)+\displaystyle \frac{{{\rm{\Theta }}}_{E}}{2}\right]\right\}\\ \,\begin{array}{c}=3{k}_{B}{N}_{enz}\left\{\displaystyle \frac{{{\rm{\Theta }}}_{E}f\left(T\right)}{g\left(T\right)\mathrm{ln}\left[{{\rm{\Theta }}}_{E}+2T\,\mathrm{ln}\,g\left(T\right)\right]}-\,\mathrm{ln}\left(\displaystyle \frac{{{\rm{\Theta }}}_{E}}{2T}\right)+\,\mathrm{ln}\,g\left(T\right)\right\}\end{array}\end{array}\end{eqnarray} \tag{ 127 }$$

Where $f\left(T\right)={e}^{{{\rm{\Theta }}}_{E}/T}+1$ and $g\left(T\right)={e}^{{{\rm{\Theta }}}_{E}/T}-1.$

According to E Gamsjäger [38], ${{\rm{\Theta }}}_{E}\cong \sqrt{\tfrac{3}{5}}{{\rm{\Theta }}}_{D}\cong 1200\,{\rm{K}},$ for a temperature of $T=309.65\,{\rm{K}},$ the value of the vibrational entropy contribution of the enzyme is given in (128):

$$\begin{eqnarray}&&\begin{array}{c}{S}_{v}^{\left(enz\right)}=0.0127706\,\eta {k}_{B}{N}_{sus}\end{array}\end{eqnarray} \tag{ 128 }$$

Entropy of rotation

A more precise treatment of the entropy of a large molecule, such as a protein or an enzyme, requires knowing its energy levels to be able to calculate its partition function. This has been considered a very complicated problem [39, 40], therefore we suggest the approach mentioned above, to treat it as a sphere of homogenously distributed mass.

The substrate, on the other hand, is by comparison a much smaller molecule. This can be modelled as if it was a nanometric ellipsoid with a homogeneously distributed mass $M,$ with its axes denoted by $a,\,b$ and $c.$ Taking its principal axes as the coordinated system, the tensor of inertia is diagonalized, and only three quantities are needed to describe it [41]:

$$\begin{eqnarray*}&&{I}_{1}=\displaystyle \frac{1}{5}M\left({b}^{2}+{c}^{2}\right),\,{I}_{2}=\displaystyle \frac{1}{5}M\left({a}^{2}+{c}^{2}\right),\,{I}_{3}=\displaystyle \frac{1}{5}M\left({a}^{2}+{b}^{2}\right)\end{eqnarray*}$$

The partition function resulting from the rotational degrees of freedom is denoted as such:

$$\begin{eqnarray*}&&{Q}_{r}={q}_{r}^{N}={\left[\displaystyle \frac{\sqrt{\pi }}{\sigma }{\left(\displaystyle \frac{8{\pi }^{2}{k}_{B}T}{{h}^{2}}\right)}^{3/2}\sqrt{{I}_{1}{I}_{2}{I}_{3}}\right]}^{N}\end{eqnarray*}$$

Which, defining $\alpha =\tfrac{\sqrt{\pi }}{\sigma }{\left(\tfrac{8{\pi }^{2}{k}_{B}}{{h}^{2}}\right)}^{3/2}\sqrt{{I}_{1}{I}_{2}{I}_{3}},$ can be compacted as seen in (129):

$$\begin{eqnarray}&&\begin{array}{c}{q}_{r}=\alpha {T}^{3/2}\end{array}\end{eqnarray} \tag{ 129 }$$

The resulting entropy of rotation is of the form (130):

$$\begin{eqnarray}&&\begin{array}{c}{S}_{r}=N{k}_{B}\,\mathrm{ln}\left(\alpha {e}^{3/2}{T}^{3/2}\right)\end{array}\end{eqnarray} \tag{ 130 }$$

The underlying difference between the contribution from the substrate and from the enzyme is contained in the values of $N$ and $\alpha .$

Entropy of rotation of the substrate

We suggest that, when the substate molecule rotates, the shape is similar to a spheroidal prolate. With this in consideration, it is possible to estimate its size based on the length of its $N-N$ bonds (1.346 _Å); therefore, its semiaxes would be:

$$\begin{eqnarray*}&&\begin{array}{l}a=1.346\times 12\,\mathring{\rm{A}} =16.152\times {10}^{-10}\,{\rm{m}}\\ b=c=1.346\times 8\,\mathring{\rm{A}} =10.768\times {10}^{-10}\,{\rm{m}}\end{array}\end{eqnarray*}$$

With its mass given as:

$$\begin{eqnarray*}&&{M}_{sus}=5.73314\times {10}^{-25}\,{\rm{kg}}\end{eqnarray*}$$

The moments of inertia are of the form:

$$\begin{eqnarray*}&&\begin{array}{l}{I}_{1}={I}_{a}=\displaystyle \frac{2}{5}M{b}^{2}=5.9858\times {10}^{-43}\,{\rm{kg}}\,{{\rm{m}}}^{2}\\ {I}_{2}={I}_{3}={I}_{b}=\displaystyle \frac{1}{5}M\left({a}^{2}+{b}^{2}\right)=4.3209\times {10}^{-43}\,{\rm{kg}}\,{{\rm{m}}}^{2}\end{array}\end{eqnarray*}$$

Since we have no information about its symmetry, we assume $\sigma =1,$ which gives the result:

$$\begin{eqnarray}&&\begin{array}{c}{S}_{r}^{\left(sus\right)}={N}_{sus}{k}_{B}\end{array}\,\mathrm{ln}\left(\alpha {e}^{3/2}{T}^{3/2}\right)\end{eqnarray} \tag{ 131 }$$

At a temperature of $T=309.65\,{\rm{K}},$ the value of the estimation of the entropy of rotation contribution by the substrate is given by (132):

$$\begin{eqnarray}&&\begin{array}{c}{S}_{r}^{\left(sus\right)}=20.8999{N}_{sus}{k}_{B}\end{array}\end{eqnarray} \tag{ 132 }$$

Entropy of rotation of the enzyme

For the estimation of the contribution to the entropy of rotation by the enzyme we consider the beta-lactamase, which will be modelled by a sphere of homogeneously distributed mass ${M}_{enz}.$ Its measurements are presented below:

$$\begin{eqnarray*}&&\begin{array}{l}{M}_{enz}=4.981616763\times {10}^{-23}\,{\rm{kg}}\,{R}_{enz}=\displaystyle \frac{1}{2}\left(5\,{\rm{nm}}\right)=2.5\times {10}^{-9{\rm{m}}}\\ I=\displaystyle \frac{2}{5}M{R}^{2}=5.1892\times {10}^{-49}\,{\rm{kg}}\,{{\rm{m}}}^{3}\end{array}\end{eqnarray*}$$

The resulting entropy of rotation for the enzyme is:

$$\begin{eqnarray}&&\begin{array}{c}{S}_{r}^{\left(enz\right)}={N}_{enz}{k}_{B}\,\mathrm{ln}\left(\alpha {e}^{3/2}{T}^{3/2}\right)\end{array}\end{eqnarray} \tag{ 133 }$$

At a temperature of $T=309.65\,{\rm{K}},$ the resulting entropy of rotation contribution by the enzyme is given by (134):

$$\begin{eqnarray}&&\begin{array}{c}{S}_{r}^{\left(enz\right)}=0.592884{N}_{enz}{k}_{B}=0.592884\eta {k}_{B}{N}_{sus}\end{array}\end{eqnarray} \tag{ 134 }$$

Entropy of equilibrium

Summing of all the contributions to the nondimensionalized entropy, (123), (124), (126), (128), (132), (134), the resulting expression is given in (135):

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}\displaystyle \frac{{S}_{eq}}{{k}_{B}{N}_{sus}}=80.3875+\eta 1.82029+0.0984549\left(1+\eta \right)+12.0203\\ \,+0.0127706\eta +20.8999+0.592884\eta \end{array}\end{array}\end{eqnarray} \tag{ 135 }$$

The entropy of fluctuation is made noticeable in the system once the interaction between enzyme and substrate begins. This moment in time will be labelled as $t=0.$

According to the Michaelis-Menten model, it is enough to follow the substrate and the enzyme-substrate complex to have a proper understanding of the system. Under such consideration, we will see that there is a decrease in entropy.

The computer simulation was carried out with the reaction rates shown in table 5 [42] below. As shown in previous sections, the simulation displays three stages of the chemical reaction. Figure 14 shows the functions fitted to the means of ${N}_{2}/{\rm{\Omega }}$ and ${N}_{3}/{\rm{\Omega }},$ where, as before, ${N}_{2}=S,$ ${N}_{3}=ES$ and ${\rm{\Omega }}={N}_{10}+{N}_{20}.$

Table 5. Transitions rates used in the simulation of the penicillin hydrolysis catalyzed by beta-lactamase.

${k}_{1}^{{\prime} }$	$41\times {10}^{6}$	$1/{\rm{mol}}$
${k}_{2}$	$2320$	$1/s$
${k}_{3}$	$3610$	$1/s$

Zoom In Zoom Out Reset image size

The method used to calculate the means was the following:

• 1.  
  The simulations were carried out $p$ number of times, with $p=1000.$
• 2.  
  Of the $p$ realizations, the one with the longest duration was selected. Its time was divided into $I$ intervals of equal width, with $I=100,\,1000,\,10000.$
• 3.  
  The datapoints of the $p$ realizations that fell in each interval were averaged.

The curves fitted serve as confirmation of an initial stage when the substrate decreases rapidly as the enzyme-substrate complex increases drastically towards the second stage. During the second stage, almost all enzymes are in the enzyme-complex state, in other words, the concentration of free enzymes fluctuates near zero; this stage is what is called the stationary state in biochemistry. The third stage is reached when substrate is depleted, therefore the concentration of enzyme-substrate complex decreases exponentially. The time duration of each stage and the functions fitted for the substrate and enzyme-substrate complex during each of these can be seen in table 6, the curves are shown in figure 14.

Table 6. Time functions that parametrize the substrate and enzyme-substrate complex concentrations in 1st, 2nd and 3rd stage.

	Duration (s)	$S$ fit	$ES$ fit
1st stage	$0\lt t\lt {t}_{1}=3.38985\times {10}^{-4}$	${\psi }_{1}\left(1\right)=0.980388-111.03901t$	${\phi }_{1}\left(t\right)=57.5649t$
2nd stage	${t}_{1}\lt t\lt {t}_{2}=0.013559$	${\psi }_{2}\left(t\right)=0.966669-70.56632t$	${\phi }_{2}\left(t\right)=0.01951$
3rd stage	${t}_{2}\lt t\lt {t}_{3}=0.0169$	${\psi }_{3}\left(t\right)=0.00983$	${\phi }_{3}\left(t\right)={e}^{38-3090t}$

Equations (32) were solved numerically, then $S\left(t\right)/N{k}_{B}$ was found using expression (36). The computer simulation allows us to find the entropy of fluctuations ${{\rm{\Xi }}}_{{\rm{jj}}}\left(0\right)$ but, even though we consider the precision to not be enough to quantitively determine ${{\rm{\Xi }}}_{jj}\left(t\right),$ we can provide an approximation using the following initial conditions:

$$\begin{eqnarray*}&&{{\rm{\Xi }}}_{11}\left(0\right)=0.377676\,{{\rm{\Xi }}}_{22}\left(0\right)=0.377676\,{{\rm{\Xi }}}_{12}\left(0\right)=0\end{eqnarray*}$$

Adding the approximation of the entropy of fluctuations to the previously obtained entropy of equilibrium, ${S}_{eq}/N{k}_{B},$ the result can be seen in figure 15. Some relevant values of the time evolution of ${S}_{eq}/N{k}_{B}$ are shown in table 7.

Zoom In Zoom Out Reset image size

Table 7. Numerical values of the entropy. ${S}_{eq}:$ entropy of equilibrium. ${\rm{\Delta }}{S}_{+}:$ entropy of fluctuation during $0\lt t\lt {t}_{3}.$ ${S}_{T}\left(t\gt 0\right):$ total entropy at the start of the reaction, $t\gt 0.$ ${S}_{T}\left({t}_{3}\right):$ Total entropy at ${t}_{3}$ (the end of the catalytic reaction). ${\rm{\Delta }}S:$ Decrease in entropy of fluctuation (${S}_{f}\left({t}_{3}\right)-{S}_{f}\left(t\gt 0\right)$).

${S}_{eq}/N{k}_{B}$	$115.83$
${\rm{\Delta }}{S}_{+}/N{k}_{B}$	$0.433$
${S}_{T}\left(t\gt 0\right)/N{k}_{B}$	$116.265$
${S}_{T}\left({t}_{3}\right)/N{k}_{B}$	$116.259$
${\rm{\Delta }}S/N{k}_{B}$	$-6.602\times {10}^{-3}\,$

The intriguing result seen in figure 15 and table 7 can be stated in two aspects:

• (1)  
  It is clear that the net entropy of the reaction is greater than it was before the start of the reaction, therefore it is a spontaneous process.
• (2)  
  The curve displayed in figure 15(a) shows a decrease in entropy, something that is only possible if there is an external energy source.

Before the reaction, the value of the total entropy is ${S}_{eq}/N{k}_{B}.$ Once the reaction starts (at $t\gt 0$), the total entropy suffers a change, increasing a quantity of ${\rm{\Delta }}{S}_{+}/N{k}_{B},$ so the total entropy at that point is ${S}_{T}\left(t\gt 0\right)/N{k}_{B}.$ As the reaction takes place, the entropy decreases until the system reaches a state of equilibrium once the catalysis has finished (at $t={t}_{3}$); the resulting decrease in entropy is ${\rm{\Delta }}S,$ so the total entropy at that point is ${S}_{T}({t}_{3})/N{k}_{B}.$ It is important to note that the magnitude of the decrement in entropy, ${\rm{\Delta }}S,$ is lower than the increase in entropy due to the reaction taking place, ${\rm{\Delta }}{S}_{+};$ in other words, ${{\rm{\Delta }}S}_{+}\gt {\rm{\Delta }}S.$ According to the second law of thermodynamics, this decrement is evidence of the existence of work during the process of catalysis.

The value of ${\rm{\Delta }}S$ is approximately 1.52% of the initial entropy of fluctuations and it is 6.71% of the mixed entropy, this could be (as stated before) because the system is receiving energy in the form of work from an external source. A revision of the existing literature leads us to suggest that it comes from the vibrational degrees of freedom of the enzyme, which will be discussed below. For processes at constant volume, the fundamental equation of thermodynamics establishes:

$$\begin{eqnarray*}&&{\rm{\Delta }}U=T{\rm{\Delta }}S+{\mu }_{2}{\rm{\Delta }}{N}_{2}+{\mu }_{3}{\rm{\Delta }}{N}_{3}\end{eqnarray*}$$

The condition of ${\rm{\Delta }}S\lt 0$ leads to the following inequality: ${\rm{\Delta }}U\lt {\mu }_{2}{\rm{\Delta }}{N}_{2}+{\mu }_{3}{\rm{\Delta }}{N}_{3}.$ Considering the quasi-stationary state as representative of the majority of the process of catalysis, then ${\rm{\Delta }}{N}_{3}=0,$ therefore the inequality takes the form of:

$$\begin{eqnarray*}&&{\rm{\Delta }}U\lt {\mu }_{2}{\rm{\Delta }}{N}_{2}\end{eqnarray*}$$

If ${\rm{\Delta }}S\lt 0$ one must have ${C}_{p}=T{\left(\tfrac{{\rm{\Delta }}S}{{\rm{\Delta }}T}\right)}_{T}\lt 0,$ so that this decrease in the entropy is consistent with results previously published: in enzyme catalysis happens that ${C}_{p}\lt 0$ [43, 44]. Also, the negative value in this heat capacity is a condition to show that in enzyme catalysis there is an optimal temperature $T={T}_{op},$ where the most efficient catalysis occurs.

The results from our computer simulation of a Michaelis-Menten system allows us to make a prediction that points towards the correct direction, since it reproduces the inequality of ${\rm{\Delta }}S\lt 0$ for the entropy of fluctuations.

However, our model has a limitation in its quantitative aspect. According to Hobbs et al [43] and Arcus et al [44], the values of the heat capacity at constant pressure are within the range $-12\tfrac{{\rm{kJ}}}{{\rm{mol}}\,{\rm{K}}}\lt {C}_{p}\lt -1\tfrac{{\rm{kJ}}}{{\rm{mol}}\,{\rm{K}}}.$ If one is to take the right-hand value of the range, at $T=298.5\,{\rm{K}},$ with ${\rm{\Delta }}T=1\,{\rm{K}},$ converting to ${\rm{eV}}$ per particle, the resulting value of the change in entropy would be:

$$\begin{eqnarray*}&&{\rm{\Delta }}S={C}_{p}\displaystyle \frac{{\rm{\Delta }}T}{T}=-3.4764\times {10}^{-5}\,{\rm{eV}}\,{{\rm{K}}}^{-1}\end{eqnarray*}$$

Furthermore, according to our results, the entropy of fluctuations per particle is:

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{\Delta }}S}{N}=-5.6875\times {10}^{-7}\,{\rm{eV}}\,{{\rm{K}}}^{-1}\end{eqnarray*}$$

This leads us to conclude that there is some aspect missing in this model.