Onsager's variational principle in active soft matter

Abstract

Onsagers variational principle (OVP) was originally proposed by Lars Onsager in 1931 [L. Onsager, Phys. Rev., 1931, 37, 405]. This fundamental principle provides a very powerful tool for formulating thermodynamically consistent models. It can also be employed to find approximate solutions, especially in the study of soft matter dynamics. In this work, OVP is extended and applied to the dynamic modeling of active soft matter such as suspensions of bacteria and aggregates of animal cells. We first extend the general formulation of OVP to active matter dynamics where active forces are included as external non-conservative forces. We then use OVP to analyze the directional motion of individual active units: a molecular motor walking on a stiff biofilament and a toy two-sphere microswimmer. Next we use OVP to formulate a diffuse-interface model for an active polar droplet on a solid substrate. In addition to the generalized hydrodynamic equations for active polar fluids in the bulk region, we have also derived thermodynamically consistent boundary conditions. Finally, we consider the dynamics of a thin active polar droplet under the lubrication approximation. We use OVP to derive a generalized thin film equation and then employ OVP as an approximation tool to find the spreading laws for the thin active polar droplet. By incorporating the activity of biological systems into OVP, we develop a general approach to construct thermodynamically consistent models for better understanding the emergent behaviors of individual animal cells and cell aggregates or tissues.

---

1 Introduction

Active matter represents a novel type of nonequilibrium systems that contain a large number of self-propelling particles or creatures moving in fluids or more complex environments.^1–10 The self-propelling units are considered to be active in the sense that they are capable of continuously converting fuel or chemical energy (stored internally or in the ambient) into directional motion or mechanical work. Active matter containing self-propelled units is ubiquitous in biology and in many artificial systems. Examples in biology are abundant and occur at all length scales, ranging from bacteria suspensions^3,9–12 to animal cells,^2,4 animal cell aggregates (or tissues),^13–17 bird and fish flocks,^18,19 and pedestrian crowds.^20,21 Artificially made active matter³ includes layers of vibrated granular rods,^1,22 collections of robots,²³ and suspensions of colloidal or nanoscale particles^24–26 that are propelled through catalytic activities at their surfaces.

A distinctive feature of active matter is that the system is locally driven out of equilibrium by active units at the length scale of a constituent component.^1–3,27 This is distinct from those nonequilibrium systems that are driven at the system boundaries.²⁸ The presence of self-propelled units in active matter breaks the detailed balance and time-reversal symmetry (TRS),^1,3,5,8 resulting in a wealth of intriguing macroscopic structures and behaviors, such as spontaneous flows,^3,4 motility-induced phase separation,^3,6,8 unusual mechanical and rheological properties,^2,10 wave propagation and sustained oscillations even in the absence of inertia,^29–31etc. One of the most interesting questions in the nonequilibrium dynamics of active matter is how the local driving forces operating at the small scale of individual active unit can produce the observable macroscopic emergent phenomena at the large scale of the whole system. Answering this question will not only shed new light on the fundamental statistical mechanics,^3,6,8,26 but also deepen our understanding of biological processes,^1–4 and help design new generations of biomimetic active materials that balance structural flexibility and stability.^1,3,5

The study of active matter can be brought into the framework of condensed matter physics based on the consideration that the collective behaviors of active matter emerge from the interactions among the constituent self-propelling units and the dissipation mechanisms operating inside the system. In particular, soft condensed matter physics provides many useful model systems for reference to active matter,^2,5e.g. the wetting of substrates by liquid droplets,^32–34 the dynamics of colloid suspensions,^35,36 the dynamics of nematic liquid crystals,³⁷ the dynamics and rheology of polymer gels,^38,39 and the phase segregation of surfactants,^40,41etc. The major challenge is to couple these model systems with active and molecularly specific processes, such as the active force generation by self-propelling units, and the binding and unbinding of transmembrane adhesion receptors on solid substrates.^2,5 Over the last decade, several soft matter systems have been revisited with a focus on this point of view.² The physical understanding for the emergent structures and behaviors of active soft matter has been rapidly growing,^1–10 with particular attention paid to dry active matter,^3,18 active polar fluids,⁴² active nematics,⁴³ active gels,⁴⁴ and active membranes.⁴⁵

Theoretically there have been two major approaches to the study of active soft matter: particle-based models^{3,6–8,26,46} and continuum phenomenological models.^{1,2,4–6,14,15,43} In particle-based models, the active units are usually modeled as self-propelled particles with fixed or variable speed and random orientation moving in an inert background, following the seminal work of Vicsek et al.⁴⁷ It provides a straightforward approach to the study of active soft matter with an emphasis on the order and fluctuations rather than the forces and mechanics.^3,7,8,26 In continuum phenomenological models, active units are represented by a smooth density field rather than individually resolved particles. A continuum model for active soft matter is usually constructed by modifying the dynamic model of a proper reference soft matter system.^1,2,4–6,43 This is typically accomplished by adding a minimal set of extra terms that cannot be derived from any free energy or dissipation functions. This is an effective way to introduce the activity and break the TRS such as active forces, active fluxes, and active chemical potentials.^1,43,48,49 However, there is another more systematic way of including activity by introducing the mechanochemical coupling between passive dissipative processes and some relevant biochemical reactions^1,4 in Onsager's framework of irreversible thermodynamics.^50,51 The two theoretical approaches are complementary. The particle-based approach involves only a small number of parameters for each active unit, and therefore the theoretical predictions can be readily compared with experiments for some model active systems such as self-propelled colloids.⁵² However, the model for interacting self-propelled particles sometimes oversimplifies the problem, and hence may lose some generality and applicability of its conclusions when applied to real systems, especially in vivo biological systems.^1,2,4,43 By contrast, the formulation of phenomenological models is based upon symmetry consideration, conservation laws of mass, momentum, and angular momentum, and laws of thermodynamics. This gives the continuum approach a large range of applicability and generality when applied to real biological processes.^1,2,4,43 In this work, we focus on the continuum phenomenological models and show that Onsagers variational principle, which has been widely used in the study of soft matter dynamics,⁵³ can be extended for the study of active soft matter.

Onsager's variational principle (OVP) was originally proposed by Lars Onsager in his seminar papers in 1931.^54,55 He showed that for irreversible processes in a near-equilibrium thermodynamic system, the thermodynamic fluxes can be written as linear combinations of conjugate thermodynamic forces, and the proportionality coefficient matrix must be positive-definite and symmetric according to Onsager's reciprocal relations (ORR). The ORR lay the foundation for the theoretical framework of linear irreversible thermodynamics. In the end of his paper, Onsager proposed OVP as a variational principle that is equivalent to the linear force–flux relations in describing dissipative dynamics. In addition, OVP can be regarded as an extension of “the principle of the least dissipation of energy” proposed by Lord Rayleigh.⁵⁶ For isothermal systems, OVP takes a simple form as follows.^53,57,58 The irreversible processes described by the thermodynamic fluxes [small alpha, Greek, dot above] follow the dynamic path that minimizes the function:

[scr R, script letter R]([small alpha, Greek, dot above]) = Φ([small alpha, Greek, dot above]) + Ḟ([small alpha, Greek, dot above];α). (1)

Here the function [scr R, script letter R] is called the Rayleighian as suggested by Doi and Edwards,⁵⁹Φ is the dissipation function which is quadratic in [small alpha, Greek, dot above] when the system is close to equilibrium, and Ḟ is the rate of change of the free energy in the isothermal system. Onsager later used his principle to study the diffusion in electrolyte solutions.⁵⁷ However, OVP has not been widely recognized and applied to describe irreversible processes for a long time until 1953 when Onsager and Machlup⁶⁰ established the statistical mechanical foundation of OVP. Since then, OVP and its relationship with other thermodynamic variational principles have been extensively studied and identified.⁶¹ In recent years, OVP has been widely used as an indispensable and powerful tool for the study of nonlinear and nonequilibrium phenomena of soft matter.^{53,58,62–65}

OVP can be used to derive many transport equations for soft matter dynamics,^53,62e.g., the Stokes equation for incompressible low-Reynolds-number flows,⁵⁸ the diffusion equation,^53,62 the reaction-diffusion equations for (low molecular weight) multi-component solutions,^51,64,65 the thin film evolution equations,^63,66–68 the phase field model for two-phase hydrodynamics,^58,69 the electrorheological hydrodynamic equation,⁷⁰ the two-fluid model for the phase separation dynamics in colloids and polymers,^53,62,63 the dynamics of polymer gels,^39,63 the dynamic equations of lipid membrane,^71–73etc. Moreover, OVP also provides a very convenient way to derive thermodynamically consistent boundary conditions that supplement the transport equations in the bulk region.^58,69 Examples include the generalized Navier boundary condition (GNBC) for the contact line hydrodynamics,⁵⁸ the generalized nemato-hydrodynamic boundary conditions for liquid crystals,⁷⁴ and the boundary conditions for block copolymer solution films,⁷⁵etc. In addition, Doi and his collaborators have recently proposed that OVP can be used as a direct variational method to find approximate solutions for complex soft matter dynamics.^63,68,76 This approximation method has been successfully used to study the evolution of droplets and thin films,^63,67,77,78 the dynamics of the beads-on-string structure of viscoelastic polymer filaments,^63,79 the sedimentation in colloidal suspensions,⁸⁰ and the translocation of a vesicle through a narrow hole,⁸¹etc. Although OVP has been widely applied with great successes in the study of inert soft matter dynamics, it is rarely used in the study of active soft matter dynamics.^65,82–84 In the present work, we will show that OVP can be readily extended to include biochemical activity and conveniently applied to study the emergent structures and behaviors of active soft matter. OVP can not only be applied to formulate thermodynamically consistent models, but also be used to generate approximate solutions for the complex dynamics of active soft matter.

This paper is organized as follows. In Section 2, a brief review of OVP is provided, and a simple extension of OVP is presented for applications to active matter, in which the active forces are treated as non-conservative forces that cannot be derived from any free energy and dissipation functions. In the next three sections, we apply OVP and its extended form to three representative active matter problems motivated by the biology of bacteria and animal cells. In Section 3, we present the first application of OVP to the directional motion of an individual active unit, e.g., a molecular motor walking on a stiff biofilament and a toy two-sphere microswimmer moving in a viscous fluid. In Section 4, we consider the two-phase hydrodynamics for active polar droplets. We use OVP to formulate a diffuse-interface model for an active polar droplet on a solid substrate. This hydrodynamic model is thermodynamically consistent and consists of hydrodynamic equations in the bulk region and boundary conditions at the solid surface. In Section 5, we consider a thin active polar droplet moving on a solid substrate in two dimensions. Under the lubrication approximation, we firstly apply OVP to derive the classical thin film equation that has been obtained previously. We then use OVP as an approximation tool to find the scaling laws for the spreading of a thin active droplet in the respective limits of negligible activity and strong activity. In Section 6, we summarize our major results, make some general remarks, and envision a few potential applications of OVP to more realistic biological problems.

2 Variational principles for active matter dynamics

In this section, we show that the original variational principle of Onsager^54,60 can be easily modified to study the dynamics of active soft matter by including active forces as non-conservative forces that can not be derived from any free energy function. Since we are mostly interested in the flow, diffusion, and biochemical reactions in active soft matter, we, therefore, limit our discussions to isothermal systems where temperature is assumed to be constant.

2.1 Onsager's variational principle and Onsager–Machlup variational principle

Consider a non-equilibrium isothermal system that is characterized by a set of coarse-grained, slow state variables, α ≡ {α₁,α₂,…,α_N}. The dynamics of the system is then described by the time evolution of α(t) which is governed, in the linear response regime near equilibrium, by the general (overdamped) Langevin equation⁶⁰in which ζ = {ζ_ij} is the friction coefficient matrix that generally depends on the state variables⁶⁵α. Here we assume that all the state variables α have the same time parity, and therefore, according to Onsager's reciprocal symmetry, the friction matrix ζ is not only positive definite (due to the second law of thermodynamics) but also symmetric,⁵⁴i.e., ζ_ij = ζ_ji. State variables with different time parities have been discussed briefly in the Appendix Section A.2. The stochastic force f_r(t) is assumed to be an uncorrelated white noise

〈f_ri(t)〉 = 0, 〈f_ri(t)f_rj(t′)〉 = 2ζ_ijk_BTδ(t − t′), (3)

with k_B denoting the Boltzmann constant and T the temperature. The generalized force f_i, in general, includes two types of forces: f_i(α,t) = f_ci(α) + f_ai(α,t), in which the conservative force f_ci(α) can be derived from some free energy F(α) byand f_ai(α,t) is the active force, a non-conservative force that cannot be derived from any energy function. Physically, the active forces arise from the persistent consumption of chemical energy and they continuously drive the system out of equilibrium locally at the small scale of individual active unit. For example, the active forces can be generated by biochemical reactions such as ATP hydrolysis in animal cells/tissues,^1,4 or by external fields such as light acting on active colloids with photosensitive coatings.^3,7 The active forces also break the time-reversal symmetry (TRS) of the system (in a sense different from the breakdown of TRS due to friction).^8,43,48,85 For example, the self-propelling force in active colloids drives the persistent motion of colloids in some direction, resulting in the intrinsic breakdown of TRS.⁴³

From the Langevin eqn (2), we calculate the transition probability⁶⁰P(α′,t + dt|α,t) from the state α at time t to α′ at t + dt:

for α′ close to α, and dt is an infinitesimal time interval. Here the Onsager–Machlup function [scr O, script letter O] is defined bywith the rates and μ_ij being the mobility coefficient matrix that is positive definite and symmetric, and satisfies ζ_ijμ_jk = δ_ik. From the transition probability P(α′,t + dt|α,t) in eqn (5), the most probable transition occurring between α and nearby α′ is the one which minimizes the Onsager–Machlup function [scr O, script letter O] or equivalently the function

[scr R, script letter R]([small alpha, Greek, dot above];α) = Φ([small alpha, Greek, dot above],[small alpha, Greek, dot above]) + Ḟ([small alpha, Greek, dot above];α) − Ẇ_a([small alpha, Greek, dot above];α) (7)

with respect to α′ or equivalently, the rates [small alpha, Greek, dot above] for prescribed α. Here [scr R, script letter R] is called the Rayleighian,^62,63 is the positive-definite dissipation-function, is the change rate of free energy, and Ẇ_a([small alpha, Greek, dot above];α) = f_ai(α,t)[small alpha, Greek, dot above]_i is the work power done by the active force, f_a. The Euler–Lagrange equation for minimizing [scr R, script letter R] in eqn (7) with respect to [small alpha, Greek, dot above] isas expected from the Langevin eqn (2) for dynamic processes where the stochastic forces and fluctuations can be neglected.† The dynamic equations in eqn (8) state the balance among the dissipative (frictional) forces, the conservative forces derived from the free energy, and the active forces that are non-conservative and generated by biochemical reactions. This variational principle of minimizing the Rayleighian [scr R, script letter R] with respect to the rates [small alpha, Greek, dot above] was originally proposed by Lars Onsager in 1931^54,60 and usually referred to as The principle of least dissipation of energy or simply Onsager's variational principle (OVP).

We would like to give some remarks on OVP as follows.

(i) A term appearing in the Onsager–Machlup function [scr O, script letter O] in eqn (6) does not contribute to the Rayleighian [scr R, script letter R] because we are considering the most probable state α′ at an immediate future time close to t, which is to be determined from the prescribed state α and hence the prescribed forces f(α,t) and Ψ(f(α,t),f(α,t)).

(ii) The active matter may also be subject to the influence of some external forces that do not arise locally from the consumption of chemical energy of the system. This can be treated by subtracting the work power Ẇ_ext([small alpha, Greek, dot above];α) done by the external forces from the Rayleighian [scr R, script letter R] in eqn (7):

[scr R, script letter R]([small alpha, Greek, dot above];α) = Φ([small alpha, Greek, dot above],[small alpha, Greek, dot above]) + Ḟ([small alpha, Greek, dot above];α) − Ẇ_a([small alpha, Greek, dot above];α) − Ẇ_ext([small alpha, Greek, dot above];α). (9)

Minimization of [scr R, script letter R] leads to an Euler–Lagrange equation that includes the external forces, in a form generalized from eqn (8).

(iii) In a continuum model of active matter, the set of slow variables α can represent the field variables both in the bulk and at the boundary. Then eqn (8) derived from OVP gives the dynamic equations both in the bulk and at the boundary, with the latter becoming dynamic boundary conditions. Furthermore, if there are external forces applied at the system boundary, then their contributions may be described by Ẇ_ext in the Rayleighian in eqn (9).

(iv) It is important to note that OVP is a local principle that can be used to find the most probable state only in the immediate future (without additional constraints). To locate the most probable paths that can take the system to the far future under various constraints,^60,87 we can divide the long time interval (say, t − t₀) into N_t sub-intervals with dt = (t − t₀)/N_t, then

where [scr O, script letter O][α(t)] is called the Onsager–Machlup functional or integral defined byor equivalently, into the following quadratic formwith being the actual rates of the system at state α and time t. Note that [small alpha, Greek, dot above] is defined by the time derivative of the state variables α, while [small alpha, Greek, dot above]* is defined by the force f(α,t) and depends only on the state variable α and time t. The Onsager–Machlup integral is non-negative definite and equals to zero (the minimum) only when [small alpha, Greek, dot above] equals to the actual kinetic path, i.e., [small alpha, Greek, dot above] = [small alpha, Greek, dot above]*. Hence the variational principle can be stated that nature chooses the kinetic path which minimizes the Onsager–Machlup functional [scr O, script letter O][α(t)] with respect to the state variables α(t). This variational principle is called Onsager–Machlup variational principle (OMVP).^63,68 In contrast to OVP, a local principle that predicts the most probable state in the immediate future, the OMVP is a global principle that can determine the most probable path taking the system to the far future under various constraints. For example, it can be used to determine the long-time behaviors such as the steady-states.^63,68 It can also be used to locate the most probable transition pathway that takes the system from one free energy minimum to another.^87,88

2.2 Direct variational methods for approximation solutions

The variational principles introduced previously can be used to obtain thermodynamically-consistent dynamic equations as well as matching boundary conditions. However, the resulting equation systems are usually difficult to solve both analytically and numerically. This subsection deals with some direct variational methods of finding approximate solutions that are obtained directly from above variational principles for the system dynamics in both short and long time scales.^63,68,76,89

Variational principles have been proposed in various fields of physics and several variational methods have been developed accordingly to find approximation solutions such as Ritz method and the least-squares method.⁸⁹ Recently, Doi developed a Ritz-type variational method based on OVP⁷⁶ by assuming that state variables α(t) is a certain function of a small number of parameters denoted by a = (a₁,a₂,…), i.e., α(t) = α(a(t)). Then the rate of the state variables [small alpha, Greek, dot above] can be written as

and the Rayleighian is then a function of the rates of parameters [scr R, script letter R](ȧ). The dynamics of the system described by the temporal evolution of the parameters a(t) is determined by minimizing [scr R, script letter R] with respect to ȧ.

Similarly, direct variational method can also be developed based on the OMVP to approximate the long time kinetic paths or states of the system.⁶⁸ We consider certain kinetic path which involves a parameter set α(a(t)). The best guess for the actual path is the path which gives the smallest value of the Onsager–Machlup functional [scr O, script letter O][α(a(t))] with respect to the parameter functions a(t). This variational method is similar to the least-square method but with a target function that is more physically meaningful based on physical principles.^68,89

These direct variational methods are useful particularly when we have an idea for the probable kinetic path and can write down the functions α(a(t)). It has been applied successfully to many problems in soft matter dynamics. In this work, we will show that these approximation methods can also be used to study the dynamics of active soft matter.

2.3 Advantages of variational approaches

The above variational principles are equivalent to Onsagers kinetic equations in eqn (8) with kinetic coefficients satisfying Onsager's reciprocal relations, but these variational principles and the direct variational methods have several advantages to investigate the dynamics of non-equilibrium systems as follows.

• Scalar formulation. The variational principles involve only physical quantities that can be defined without reference to a particular set of generalized coordinates, namely the dissipation function, free energy, and active work power. This formulation is therefore automatically invariant with respect to the choice of coordinates for the system, which allows us a great flexibility in choosing state variables and rates.

• Thermodynamic consistency. The variational principles incorporate the intrinsic structure of Onsagers theory of non-equilibrium thermodynamics clearly. They provide compact invariant ways of obtaining thermodynamically-consistent dynamic equation systems where the pairs of rates and forces are obtained automatically.

• Direct variational approximation tools. The direct Ritz-type variational method of finding approximation solutions for the system dynamics bypasses the derivation of the Euler–Lagrange equations and goes directly from a variational statement of the problem to the solution of the Euler–Lagrange equations. This approximation method helps to pick up the most important dynamic behaviors and to simplify the calculations significantly from complicated partial differential equation systems to simple ordinary differential equations. In addition, the direct least-square-type variational method based on OMVP of minimizing the Onsager–Machlup integral further optimizes the search for more realistic kinetic paths and provides a new method of studying long-time steady-state dynamics of the system.

The above variational principles have been successfully applied to the diffusion in electrolyte solutions by Onsager himself in 1940s,⁵⁷ and more recently applied to various soft matter systems such as multiphase flows,^58,69 electrorheological fluids,⁷⁰ colloid suspensions,⁸⁰ polymer solutions,^53,62,63 polymer gels,^39,63 liquid crystals,^62,74 vesicles,⁸¹ membranes^71–73 and so on. This indicates that OVP is an important principle in soft matter dynamics.^62,63,68 In this work, we present its applications to active soft matter dynamics that is mostly motivated by biological applications.

Before ending this section, we would like to summarize the general steps for applying OVP to the dynamics of active soft matter for which the dynamic equations are not yet known or still controversial.^62,63

(i) Choose a set of coarse-grained, slow variables, α ≡ {α₁,α₂,…}, to describe the time evolution of the macroscopic state of the system.

(ii) Construct the free energy function, F(α), and calculate the rate of change of the free energy, Ḟ([small alpha, Greek, dot above];α).

(iii) Construct the dissipation function, Φ([small alpha, Greek, dot above],[small alpha, Greek, dot above]), which is quadratic in the rates/fluxes [small alpha, Greek, dot above].

(iv) Find the work power done by the active forces, Ẇ_a([small alpha, Greek, dot above];α), based on the specific activity considered, and find the work power done by some other external forces, Ẇ_ext([small alpha, Greek, dot above];α). The external forces are usually applied at the system boundary and do not arise locally from the consumption of chemical energy of the system.

(v) Minimize the Rayleighian in eqn (9): [scr R, script letter R]([small alpha, Greek, dot above];α) = Φ([small alpha, Greek, dot above],[small alpha, Greek, dot above]) + Ḟ([small alpha, Greek, dot above];α) − Ẇ_a([small alpha, Greek, dot above];α) − Ẇ_ext([small alpha, Greek, dot above];α), with respect to the rates/fluxes [small alpha, Greek, dot above]. Note that some additional constraints on the system dynamics may need to be imposed by using Lagrange multipliers.

Furthermore, if we have an idea about the most probable kinetic path, then we can write down the slow variables α = α(a(t)) as functions of a small number of parameters, a = (a₁,a₂,…). We can follow the above steps and obtain the Rayleighian as a function of ȧ and a as [scr R, script letter R](ȧ;a). The minimization of [scr R, script letter R] with respect to ȧ will then provide an approximate description for the active matter dynamics directly.

3 Applications 1: directional motion of individual active units

Activity in biology or in some artificial active systems usually arises from the consumption of ATP and the mechanochemical cross-coupling, but sometimes activity appears simply as time-dependent constraints in geometric shapes. In this section, we consider the first applications of OVP to the directional motion of an individual active unit as shown in Fig. 1: a molecular motor walking on a polar biofilament and a toy two-sphere microswimmer swimming in a viscous fluid.

Fig. 1 Directional motion of individual active units. (a) A myosin motor catalyzes ATP hydrolysis and converts the released chemical energy into its directional motion on an actin filament⁹⁴ toward the plus (barbed) end. Meanwhile, an external force, f_ext, is applied on the myosin through optical trapping of a nano-bead attached to the motor. (b) A toy two-sphere microswimmer swims in viscous fluids. The body length, 2[small script l](t), of the microswimmer oscillates cyclically and the friction coefficient, ζ, is asymmetric for forward motion (with smaller friction) and backward motion (with larger friction), as shown in eqn (23). The front-back asymmetry in friction is indicated by inclined thorns on the microspheres.

3.1 Molecular motors walking on biofilaments: mechano-chemical cross-coupling

Molecular motor proteins are enzymes that bind adenosine triphosphate (ATP) and catalyze its hydrolysis to adenosine diphosphate (ADP) and inorganic phosphate (P_i):⁹⁰

ATP ⇌ ADP + P_i. (14)

The chemical energy released from this ATP hydrolysis is partially converted into mechanical work or directional motion of motors along some stiff biofilaments that are made of other proteins. Animal cells in vivo contain over a hundred different motor proteins, which can be classified into three different families: myosins moving along actin filaments, kinesins and dyneins moving along tubulin filaments. The underlying biofilaments are usually periodic and fairly rigid structures with a period ξ ∼ 10 nm. They are moreover polar or asymmetric, so that one can define a “plus” end and a “minus” end. A given motor always moves in the same deterministic direction: myosin moves along actin filaments towards their plus end, and kinesins and dyneins move along tubulin filaments towards their plus and minus ends, respectively. Motor molecules play a key role in cell contraction, cell division, intracellular transport, and material transport along the axons of nerve cells, etc.⁹⁰

In this subsection, we use OVP to formulate a thermodynamic description in the linear (near-equilibrium) regime for the directional motion of a translationary molecular motor along a polar filament against an external force, as shown in Fig. 1a. This Onsager-type description is pioneered by Kedem & Caplan⁹¹ and extended by Chen & Hill.⁹² We take our thermodynamic system to include the molecular motor and the surrounding solution of ATP, ADP, and P_i. The system is coupled to a heat reservoir and a work reservoir, which can apply external forces on the motor, for example, by optical tweezers, that is, optical trapping of a nano-probe attached to the motor (see Fig. 1a). The states of the thermodynamic system can be described by the average motor position x, the polarization vector p (describing the polarity of the filament and assuming to point from minus end to plus end), and the average number N_α of chemical components involved in ATP hydrolysis with α = ATP, ADP, and P_i.

The reaction free energy for ATP hydrolysis takes the form of

[scr F, script letter F]_r = [scr F, script letter F]_r(N_ATP,N_ADP,N_Pi), (15)

from which we find the rate of the change of free energy as⁴Here r ≡ dξ/dt is the reaction rate of ATP hydrolysis with ξ being the reaction extent, dξ ≡ dN_α/ν_α = −dN_ATP = dN_ADP = dN_Pi, and ν_α being the stoichiometric coefficients (negative for reactants and positive for products, here ν_ATP = −1, ν_ADP = +1, and ν_Pi = +1). Δμ is the reaction affinity of ATP hydrolysis, given bywhich measures the free-energy change for the hydrolysis of each ATP molecule. At chemical equilibrium Δμ = 0, whereas it is positive when ATP is in excess and negative when ADP is in excess. Under in vivo conditions, ATP is usually in excess with Δμ > 0; the reason may be that ATP has evolved as a biological hydrotrope to keep biomolecules soluble at high concentrations and subsequently used a the “energy currency” of the cell due to its high energy phospho-diester bonds.⁹³

The irreversible dynamics of the thermodynamic motor/filament system is characterized by two rates: the reaction rate r and the average motor velocity v = ẋ. In the linear response regime close to equilibrium, the dissipation function is a quadratic function of the rates given by

and the rate of work done by the external force f_ext to the motor is given by

Ẇ_ext = f_ext·v. (19)

Here Λ and ζ are generalized friction coefficients, and λ is the mechanochemical coupling coefficient, which is nonzero only if the filaments are polar and can be either positive or negative.⁹⁵ The positive-definiteness of the dissipation function requires Λ,ζ > 0 and Λζ − λ² > 0. The degree of mechanochemical cross-coupling can be quantified by (with −1 < q < 1), as suggested by Kedem and Caplan.⁹¹

Using the Rayleighian in the presence of the external force and eqn (16)–(19), we minimize [scr R, script letter R] with respect to r and v and obtain

Δμ = Λr − λv, (20a)

f_ext = −λr + ζv. (20b)

Here we have taken v and f_ext to be in parallel with p; a positive v (or f_ext) means the motion of the motor (or the direction of the external force) is along the direction of p, pointing toward plus end as assumed. Note that the Onsager reciprocal relation for mechanochemical cross-coupling is automatically satisfied. Physically, eqn (20b) can be written as a balance equation f_ext + f_a + f_fric = 0 for the external force f_ext, the active force f_a = λr, and the frictional force f_fric = −ζv.

To be specific, for the in vivo motion of myosin motors along actin filaments, ATP is usually in excess with a constant Δμ > 0 and myosin motor always move towards to the plus end of the actin filament. Therefore, the load-free motor velocity (for f_ext = 0) must be positive (and hence λ > 0 and 0 < q < 1). In this case, we can identify the following four regimes from eqn (20):

(i) For f_ext > 0, we have v > 0 and r > 0 (hence f_extv > 0 and rΔμ > 0), the external force pulls the motor to the plus end. Meanwhile, the excess ATP hydrolyzes and the chemical energy is consumed to drive the motor along the same direction to the plus end.

(ii) For small negative force f_ext < 0 and f_ext > f_stall,v, we still have v > 0 and r > 0 (hence f_extv < 0 and rΔμ > 0), the excess ATP hydrolyzes and the released chemical energy is converted into mechanical work. Here f_stall,v = −λΔμ/Λ < 0 is called stall force of myosin motion and when f_ext = f_stall,v, the motor is stationary with v = 0.

(iii) For f_stall,r < f_ext < f_stall,v, the moving direction of the motor is reversed with v < 0 and hence f_extv > 0. That is, the external force is doing positive work on the motor moving along the plus end. However, we still have r > 0 and hence rΔμ > 0. That is, the excess ATP hydrolyzes and the released chemical energy also drives the motion the motor towards the same plus end. Here f_stall,r = −Δμζ/λ < 0 is called stall force of ATP hydrolysis and when f_ext = f_stall,r, the ATP hydrolysis is inhibited with r = 0.

(iv) For f_ext < f_stall,r < 0, we have v < 0 and r < 0, hence f_extv > 0 and rΔμ < 0. That is, the external force is doing positive work on the motor and produces ATP that is already in excess; the system then works as an ATP pump.

Therefore, the motor/filament system is a reversible machine: it can not only convert chemical energy into mechanical work, but can also convert mechanical work into chemical energy. In this work, we are particularly interested in the regime (ii), in which ATP hydrolysis occurs spontaneously and the released chemical energy is used to drive the system out of equilibrium continuously.

Now let's consider a practical limit at which the ATP hydrolysis rate r (or equivalently the active force f_a = λr) is taken as a given positive parameter that measures the activity of the system. That is, the effect of mechanical forces on the ATP hydrolysis is neglected and the rate of ATP hydrolysis is determined dominantly by the chemical affinity as r ≈ Λ⁻¹Δμ > 0. This leads to a reduced description in which the position of the molecular motor in directional motion becomes the only state variable, while the amounts of the reactants and products in the ATP hydrolysis are no longer involved. As a result, the Rayleighian reduces to its extended form in eqn (9) as

in which is the dissipation function, f_av = Ẇ_a is the rate of work done by the active force to the reduced system, and f_extv = Ẇ_ext is the rate of work done by the external force. Minimizing this Rayleighian gives the force balance equation −ζv + f_a + f_ext = 0.

Finally, we would like to point out that the results obtained here from thermodynamic description for the motion of molecular motors on polar filaments are completely independent of any underlying microscopic mechanisms. However, the above linear-response theory applies only to the linear regime near equilibrium where Δμ/k_BT ≪ 1 and f_extξ/k_BT ≪ 1 with ξ being the typical molecular size of relevant proteins. In real life, molecular motors mostly operate far from equilibrium (with Δμ ∼ 10k_BT) and the velocity v(f_ext,Δμ) and the rate of ATP consumption r(f_ext,Δμ) are in general highly nonlinear. Therefore, more specific models such as a minimal two-states model for molecular motors should be constructed to arrive at a more comprehensive understanding of the specificity and robustness of the directional motion of motors in highly fluctuating environment.⁹⁵

3.2 A toy two-sphere microswimmer: active shape changes

Many animals and cells can actively change their shape in some periodic or cyclic manner to migrate on frictional substrates or swim in their surrounding viscous fluid environment.^9,10,96–98 For example, snakes and some worms migrate on the ground by generating body waves to change their shapes.⁹⁶ Some bacteria such as Escherichia coli, swims in fluids through bundling and rotating their flagella as driven by rotary motors.^9,10 Many types of animal cells can also migrate on substrates by dramatic periodical shape changes⁹⁷ as primarily driven by their active cellular cortex that consumes chemical energy. It is interesting to note that in most of the cell migration driven by active shape changes, the cell migration velocity shows highly nonlinear dependence on the active force or the active shape-changing velocity of the cell. This seems to go far beyond the linear-response regime and be out of the scope of OVP. However, we will show that OVP developed in the linear-response regime can still be employed if we expand the set of state variables⁶³ to include not only the center-of-mass position of the cell but also the fast changing body length which describes the shape change.

Specifically, to show how periodic shape changes can generate directional self-propulsion, here we consider a toy microswimmer that is composed of two microspheres⁹⁹ as shown in Fig. 1b. Let x₁ and x₂ denote the coordinates of these two microspheres. Then the directional motion and the shape changes of the microswimmer can be described by the temporal evolution of the center-of-mass position and the half-body-length of the swimmer , which are taken as the two slow variables. The toy microswimmer can actively change its shape by periodically changing its body length (or the distance between the two spheres) 2[small script l](t) as

[small script l](t) = [small script l]₀ + a sin(ωt), (22)

where [small script l]₀ is the half of the average body length of the microswimmer, a and ω are the amplitude and frequency of the shape oscillations, respectively, and the shape-oscillation period is T_{[small script l]} = 2π/ω.

The toy microswimmer subjected to the periodic shape oscillations can achieve directional motion only when there exists some mechanisms that break the front-back symmetry. Here we consider an asymmetry in the viscous friction, defined by f_v = −ζ(ẋ)ẋ, in which the friction coefficient ζ(ẋ) depends on the moving direction of each microsphere according to:

The dynamics of the microswimmer is characterized by the velocities (the rates), ẋ₁ and ẋ₂, of the two spheres. To the leading order in the two rates, the dissipation function is given bywhere v = ẋ_c is the center-of-mass velocity of the microswimmer, ζ₁ and ζ₂ are the frictional coefficients given in eqn (23) that depend on the signs of and , respectively. Note that the dissipation function Φ_t is actually not quadratic in dissipative rates any more but highly nonlinear; a highly nonlinear dissipation function has been discussed in classical mechanics of particles before.¹⁰⁰

In most microswimmers, it is natural to assume that there is a clear separation of time scales between their shape oscillations and the directional motion. The directional motion of the microswimmer is usually much slower than its shape oscillations, i.e., [small script l]₀/v ≫ T_{[small script l]}. We can, therefore, integrate out the relatively fast varying variable, the half-body-length [small script l](t), in one cycle of shape oscillation and arrive at a time-averaged dissipation function of the slow variable, the center-of-mass velocity v, by as

where for |v| ≤ aω and [scr A, script letter A] = 1 − ζ₊/ζ₋ is a dimensionless parameter measuring the degree of the front-back asymmetry in the friction coefficients. Note that this dissipation function is not quadratic but highly nonlinear in the rate v. We can define the Rayleighian [scr R, script letter R](v) = Φ − Φ(v = 0) to determine the time-averaged directional motion of the toy microswimmer. Minimization of [scr R, script letter R](v) gives the cell migration velocity v as a function of the active velocity aω of the periodic shape change, and this function can be compared to experimental observations. Furthermore, in the limit of small [scr A, script letter A], we have |v| ≪ aω and the Rayleighian can be approximated, to the leading order, bywhere ζ_eff = ζ₋ + ζ₊ is the effective frictional coefficient and is a time-averaged effective active force that drives the directional self-propulsion of the microswimmer. Minimization of [scr R, script letter R] gives the directional velocity of the microswimmer as

v = f_a,eff/ζ_eff, (27)

which approaches v = aω[scr A, script letter A]/π in the limit of weak asymmetry with ζ₊/ζ₋ → 1⁻ or [scr A, script letter A] → 0⁺.

We would like to give some remarks on the directional motion of the toy two-sphere microswimmer as follows.

(i) Most dynamic behaviors of biological systems show strong nonlinearity. For example, in the toy two-sphere microswimmer, the swimmer migration velocity shows highly nonlinear dependence on the active force or the active shape-changing velocity of the swimmer. However, in many cases, OVP can still be employed if we expand the set of slow state variables properly.⁶³ For example, here our set of slow variables includes not only the center-of-mass position of the swimmer but also its fast-changing body length.

(ii) The dissipation function is non-zero even for symmetric microswimmers (with ζ = ζ₋ = ζ₊) when there is no average directional motion (i.e., v = 0): . It arises in the symmetric microswimmer from the viscous dissipation due to the fast shape-oscillation in viscous fluids.

(iii) Similar to the walk of molecular motors in the previous example, the active shape changes of microswimmers are also driven by spontaneous ATP hydrolysis. Then the irreversible dynamics of the microswimmer should be characterized by the rate of ATP hydrolysis r in addition to the sphere velocities, ẋ₁ and ẋ₂. The rate of the change of free energy is given in eqn (16) by . The dissipation function is given, to the leading order in the rates, by

Also as mentioned in the previous example, in many practical cases we can take a limit at which the ATP hydrolysis rate r and hence the active force f_a = λr is a given positive parameter. This leads to a reduced description in which the dynamics of the microswimmer is described only by the sphere velocities. In this case, the Rayleighian takes to the general form of eqn (7) asin which the first two terms compose the dissipation function Φ(ẋ₁,ẋ₂), the last two terms compose the rate of work done by the active forces to the system with Ẇ_a = f_aẋ₁+ (−f_a)ẋ₂. Minimizing this Rayleighian with respect to ẋ₁ and ẋ₂ gives the force balance equations for each sphere:

f_a − ζ₁ẋ₁ = 0, −f_a − ζ₂ẋ₂ = 0, (30)

respectively. Here the pair of active forces f_a and −f_a forms a force dipole on the microswimmer. If the rate of the ATP hydrolysis r is oscillating and the resulted active force f_a = λr takes the form ofthen it will drive a shape oscillation defined in eqn (22). Note that the effective active force f_a,eff in eqn (26) is not a simple time-average of the oscillating active force f_a in eqn (31), but is the “net” active force that drives the directional motion of the microswimmer.

Finally, we would like to point out that in our toy two-sphere microswimmer, the two necessary conditions for a steady directional motion are the active shape oscillations as the energy input and the frictional asymmetry that breaks the front-back symmetry. Similarly, for a long thin swimming micro-filament, the hydrodynamic friction is anisotropic: it experiences less friction when moving along its axis than perpendicular to it. In this case, a cyclic beat pattern on the filament will be able to drive directional propulsion in a similar manner as in the above one-dimensional toy microswimmer‡.

4 Applications 2: active polar fluid models for collectives of active units

Let's now consider an active fluid that includes collectives of active units with anisotropic shapes and polarity such as rod-like self-propelled colloids,³ rod-like bacteria,^3,9,10 and the active networks of stiff filaments in the cytoskeleton of living cells.^1,2 The constituting active units can assume chemical energy to apply (extensile or contractile) dipole forces to their surrounding inert environment that drive the system locally out of equilibrium, as schematically shown in Fig. 2. Such active fluids are in contrast to more familiar passive (inert) non-equilibrium systems which are usually driven externally at their boundaries.²⁸

Fig. 2 Active stresses generated by active units: (a) extensile force dipoles generated by the bacterial microswimmer, and (b) contractile force dipoles generated by the actomyosin filament.

4.1 Active polar fluid regarded as a reactive fluid

Active fluids with polar constituent agents often show phase separation (with coexisting dilute and dense phases) and collective orientational or polar order (with collective alignment on average).^1,43 Such an active polar fluid can be viewed as a reactive fluid where viscous flows and diffusion are closely coupled with biochemical reactions.^1,4 A generalized hydrodynamic theory has been developed to describe the dynamics of such active polar fluids, phenomenologically based on conservation laws and symmetry considerations,^1,4 which is in complete compatibility with micro-reversibility and Onsager's reciprocal relations. This is in contrast to other methods of modeling active polar fluids where some active (non-equilibrium) terms are selectively added to the dynamical equations for their passive counterparts to break TRS.^43,48,85,101 In these methods, the rates of biochemical reactions in the active fluids are implicitly assumed to be a constant that is independent of the surrounding mechanical environment and simply determined by some preset reaction affinity. In this case, the active terms arising from biochemical reactions can not be derived from any free energy or dissipation function, and hence can only be added in an ad hoc manner.

In this subsection, we present an alternative derivation of the generalized hydrodynamic model for an active polar fluid that is regarded as a reactive fluid involving the ATP hydrolysis/synthesis. In the next subsection, we will show that in the same active polar fluid, if the effects of polarization and flow on the ATP hydrolysis are negligible, then the rate of ATP hydrolysis becomes a constant simply determined by the preset constant chemical affinity Δμ. The activities, driven by the spontaneous ATP hydrolysis, are then represented by local external non-conservative force fields that are added as the active terms to the dynamic model of a passive polar fluid.

To be specific, here we use OVP to derive a diffuse-interface model for a droplet of active polar fluids moving on a solid substrate, as schematically shown in Fig. 3a. The states of such an active polar droplet can be described by the following slow field variables: the scalar composition field ϕ(r,t) (distinguishing the coexisting passive isotropic phase from the active polar phase), the polarization vector field p(r,t) (describing the average orientation of active polar agents), the average fluid velocity field v(r,t), and the density field n_α(r,t) of chemical components involved in ATP hydrolysis (eqn (14)) with α = ATP, ADP, and P_i. For an active polar fluid that is confined between solid substrates or flows at the solid surfaces, the total free energy includes four contributions, [scr F, script letter F][ϕ,p,n_α] = [scr F, script letter F]_ϕ + [scr F, script letter F]_p + [scr F, script letter F]_r + [scr F, script letter F]_s, as respectively given by

Here [scr F, script letter F]_ϕ is the free energy for the two (isotropic and polar) phase coexistence, in which a and the stiffness parameter K_ϕ are both constants. [scr F, script letter F]_p is the free energy for polar liquid crystallinity, in which for simplicity, we employ the approximation of one elastic constant^37,102K_p and we take b₁(ϕ) = −b₀ϕ with b₀ > 0 that controls the isotropic-to-polar phase transition: b₁ < 0 in the polar phase (with ϕ = 1) and b₁ > 0 in the isotropic phase (with ϕ = −1). Note that the last cross-coupling term in [scr F, script letter F]_p represents the cases of perpendicular anchoring at the isotropic-polar phase interfaces and defines the orientation of p: if c > 0, p points from polar phase to isotropic phase. [scr F, script letter F]_r is the reaction free energy for ATP hydrolysis: if df_r/dn_ATP > 0 (that is, ATP is in excess), the forward ATP hydrolysis is exergonic, occurs spontaneously, and can be used to drive the changes in the motor configurations and generate mechanical motion, resulting in active dipole forces on the the surrounding passive polar fluids (see Fig. 2). In this case, the polar fluid will never reach thermodynamic equilibrium states and will be driven locally out of equilibrium by the active units or the motors that consumes ATP. [scr F, script letter F]_s is the surface energy at the substrate surfaces (with a_s and b_s being constants, and A being the surface area), which characterizes the adhesion strength of active units to the surface and the anchoring conditions for the agent orientation. Note that such a free energy [scr F, script letter F] will stabilize a droplet of active polar phase (ϕ ≃ 1 and |p| ≃ 1) in coexistence with surrounding fluids of passive isotropic phase (ϕ ≃ −1 and |p| ≃ 0), or vice versa.⁴³

Fig. 3 A droplet containing a large number of active units on a solid substrate. (a) A general two-phase diffuse-interface model of an active droplet moving on solid substrates. A droplet of active polar phase is stabilized to coexist with surrounding fluids of passive isotropic phase. (b) Spreading of a thin active droplet on a solid substrate. Planar or homogeneous anchoring condition is assumed at all bounding surfaces with which the active units are in contact.

The composition variable, ϕ, is a conserved phase parameter and its dynamics follows the following conservation equation

∂_tϕ = −∇·(ϕv + J). (33)

However, the polarization vector p is not conserved and its rate of change is defined by ṗ ≡ ∂p/∂t + v·∇p. Furthermore, the density fields n_α (with α = ATP, ADP, and P) are also not conserved due to the presence of chemical reaction of ATP hydrolysis or synthesis, and n_α follows the dynamic equation of the form

∂_tn_α = −∇·(n_αv) + rν_α. (34)

Here ν_α is the stoichiometric coefficients (see the discussion about ATP hydrolysis after eqn (16)). The reaction rate r ≡ dξ/dt depends on the concentrations of all chemical components. (This dependence is one of the constitutive relations to be derived later from OVP and given in eqn (45c).) However, in the present work, we do not intend to go into the specific expression of this concentration dependence because our purpose is to show that accompanying the chemical reaction, a mechanical force arises from the mechanochemical coupling as the active force. It is also noted that the general dynamics of ATP hydrolysis should be described by reaction-diffusion equations. However, for simplicity, here we neglect the diffusion processes, assuming that the density of each component is simply advected by the flow and produced or consumed by the chemical reaction.

Using eqn (33) and (34) and the definition of ṗ, we obtain the change rate of free energy from eqn (32) as

where σ^e is the Ericksen stress tensor^37,62,74 given by

σ^e = −p̂I − K_ϕ∇ϕ∇ϕ − K_p∇p_k∇p_k − cp∇ϕ, (36)

and σ^e satisfies the generalized Gibbs–Duhem relationHere is the generalized pressure with f̂_m (m = ϕ, p, r) being the volume density of free energy in eqn (32), μα ≡ ∂f_r/∂n_α being the chemical potential, and the reaction affinity, Δμ, of ATP hydrolysis given in eqn (17). The chemical potentials μ_ϕ in the bulk and L at the solid surfaces are, respectively, given by

L ≡ a_s + K_ϕn̂·∇ϕ + cn̂·p, (39)

with denoting the derivative of b₁. The molecular field h in the bulk and H at the solid surfaces are, respectively, given by

H ≡ b_sp·n̂n̂ + K_pn̂·∇p. (41)

The energy dissipation function Φ is a quadratic function of three dissipative rates: the shear rate , the rate of change of polarization ṗ, and the rate of ATP hydrolysis r. These rates have the same time parity and from symmetry considerations, Φ can be written into the following invariant scalar form as⁶²

where the frictional coefficients β₁,…,β₁₁ are constants and the resulted frictional coefficient matrix, β, is positive definite. That is, the diagonal coefficient terms in β are all positive, but the off-diagonal coefficient terms that describe the cross-coupling effects between two dissipative processes can be negative (although some other inequality relations have to be satisfied to keep the positive-definiteness of β). Note that the proper dissipative rates that are present in Φ and associated with the orientational dynamics of p is not ṗ but the convected co-rotational time derivative^37,62,63,74 of p:

P ≡ ṗ − ω × p = ṗ + Ω·p, (43)

which characterizes the rotation of p relative to the rotation of surrounding fluids with being the vorticity, and being the antisymmetric part of the velocity gradient tensor. Moreover, for simplicity, we have neglected the cross coupling among the transport of ϕ, the dynamics of polarization, and the ATP hydrolysis in the bulk fluid. This coupling can represent the treadmilling dynamics of the constituting components.¹⁰³ We have also ignored the possible dissipative relaxation processes associated with the anchoring of p at the solid surfaces. In addition, we would like to point out that the choice of dissipative rates in active polar fluids is not unique.^50,51 Another set of rates has been taken by Marchetti et al.:¹ the viscous stress σ^v (or the momentum flux), r and P. In this case, the time parity of σ^v is different from the other two rates r and P. A brief discussion about the consequences of different choices of dissipative rates or thermodynamic fluxes on the symmetry of Onsager coupling matrix and the applications of OVP has been presented in the Appendix Section A.2.

Then the Rayleighian is given by

where the local incompressibility constraints, ∇·v = 0, have been taken into account with the pressure p being the Lagrange multiplier. Minimization of [scr R, script letter R][v,J,ṗ,r] with respect to the rates gives the dynamic equations for two-phase active polar flows as

−∇p + ∇·(σ^e + σ^v + σ^a) = 0, (45a)

h = γ₁P + γ₂p·[small epsi, Greek, dot above] − h^a, (45b)

Δμ = −λp·P − [small zeta, Greek, tilde]pp:[small epsi, Greek, dot above] + Λr, (45c)

together with the incompressibility condition ∇·v = 0, the dynamic equations for ϕ in eqn (33) and for n_α in eqn (34). The stress tensors, σ^e, σ^v, and σ^a, are, respectively, given by eqn (36), and

σ^v = α₁([small epsi, Greek, dot above]:pp)pp + α₂pP + α₃Pp + α₄[small epsi, Greek, dot above] + α₅pp·[small epsi, Greek, dot above] + α₆[small epsi, Greek, dot above]·pp, (46a)

σ^a = −[small zeta, Greek, tilde]rpp, (46b)

with α_i being the Leslie viscosity coefficients. Note that p̂ in σ^e in eqn (36) can be absorbed into p in eqn (45a) due to the incompressibility of the active polar fluids. The active molecular field h^a is driven by ATP hydrolysis and is defined via

h^a ≡ λrp. (47)

The diffusion flux in eqn (33) is given by J = −M∇μ_ϕ. Note that the mechanochemical cross-coupling indicated in eqn (45) is similar (although more complex) to eqn (20) for the mechanochemical coupling in the motor/filament system in Section 3.1. The general discussions there about the reversible conversion of chemical energy and mechanical work also apply here in active polar fluids. We will also focus on the regime where ATP is in excess and its hydrolysis occurs spontaneously and the released chemical energy is used by the suspending active units to drive the surrounding inert polar fluids out of equilibrium continuously.

Furthermore, from the minimization of [scr R, script letter R][v,J,ṗ,r], we can also obtain the thermodynamically-consistent boundary conditions at the solid surfaces that supplements the dynamic equations in the bulk fluids:

∂_tϕ + v_τ·∇_τϕ = −ΓL(ϕ,p), (48a)

together with the impermeability boundary conditions,

n̂·v = 0, n̂·J = 0, (48c)

and the equilibrium anchoring boundary condition,

H ≡ b_sp·n̂n̂ + K_pn̂·∇p = 0. (48d)

Here the subscript τ denotes the tangential component in the plane of the substrate. Note that the boundary condition for the tangential velocity in eqn (48a) is similar to the generalized Navier boundary condition (GNBC) formulated for immiscible two-phase flows at solid surfaces.⁵⁸ The frictional coefficient β now depends on the adhesion strength of the active fluids on the substrate surfaces.

The phenomenological coefficients in the above equation systems are functions of the frictional coefficients β_i in the dissipation function Φ in eqn (42) as⁶²

with γ₁ = α₃ − α₂ and γ₂ = α₆ − α₅. Note that the Parodi relation α₂ + α₃ = α₆ − α₅ and the symmetry of coefficients for cross-coupling effects are all automatically satisfied.

4.2 Active polar fluid regarded as a passive polar fluid under local external non-conservative fields

In the last part of Section 3.1, we have considered a practical and useful limit at which the effect of mechanical forces on the ATP hydrolysis can be neglected and the rate of ATP hydrolysis is determined dominantly by the chemical affinity as a constant parameter. Here we consider the same limit of the above general dynamics of active polar fluids: the effects of polarization and flow on the ATP hydrolysis are negligibly small, and the rate r in eqn (45c) is constant and simply determined by the reaction affinity or the change of reaction free energy, i.e., r ≈ Λ⁻¹Δμ. In this case, the rate of change of the free energy and the dissipation function reduce toIt is interesting to note that the two mechanochemical cross-coupling terms (the two terms with coefficients β₇ and β₈) in [capital Phi, Greek, tilde], parametrized by the constant reaction rate r, can be rewritten as ,in which the active stress σ^a defined in eqn (46b) and the active molecular field h^a defined in eqn (47) now become

σ^a = −ζpp, (53a)

h^a = ζ_hp, (53b)

with ζ ≡ [small zeta, Greek, tilde]Λ⁻¹Δμ and ζ_h ≡ λΛ⁻¹Δμ. The cross-coupling coefficients β₇ and β₈ (hence ζ and ζ_h) can be either positive or negative as long as they can preserve the positive definiteness of the dissipation function. For example, negative and positive values of ζ correspond to contractile and extensile active stresses, respectively, as schematically shown in Fig. 2.

We would like to further point out that eqn (52) indicates that in the limit of constant reaction rate of ATP hydrolysis, the ATP-induced activity in the active polar fluid can be regarded simply as some local non-conservative fields applied externally on the passive polar fluid. The active characteristic of these external fields is reflected in the fact that the active stress σ^a and the active molecular field h^a both depend on the local state variable (the polarization), p. Furthermore, these active fields driven by spontaneous ATP hydrolysis break the time-reversal symmetry of the polar fluids.

It follows that according to the general form of eqn (7), the Rayleighian in eqn (44) can be rewritten as

Here represents the work power done on the passive polar fluids by the active stress σ^a and the active molecular field h^a. The dissipation function [capital Phi, Greek, circumflex] for the passive polar fluid is given by

Minimization of [scr R, script letter R] gives the following simplified dynamic equation system:

(i) The dynamic equations for v: the incompressibility condition ∇·v = 0, the generalized Stokes' equation in eqn (45a) with the stress tensors σ^e, σ^v, and σ^a, given in eqn (36), (46a), and (53a), respectively;

(ii) The dynamic equation for ϕ: the conservation equation in eqn (33) with J = −M∇μ_ϕ;

(iii) The dynamic equations for p:

P = γ₁⁻¹h − γ₁⁻¹γ₂p·[small epsi, Greek, dot above] + γ₁⁻¹h^a, (56)

with the active molecular field h^a given in eqn (53b) and the passive molecular field h given in eqn (40);

(iv) The boundary conditions in eqn (48) still apply to the present case.

We would like to point out that in Section 4.1, a complete model is constructed to incorporate the chemical reaction and explicitly describe the mechanochemical coupling. In this description, the time-reversal symmetry (TRS) is preserved, and so is Onsagers reciprocal relation (ORR) for mechanochemical coupling. In Section 4.2, the limit of constant reaction rate is taken, and a simplified model is obtained from the complete one in Section 4.1. In this limit, the TRS is lost, and so is ORR for mechanochemical coupling. However, the Parodi relation for the cross coupling in the passive polar fluid is still preserved.

Finally, we note that the dynamic equation system for two-phase active polar flows on solid substrates is similar to that for two-phase passive polar flows on solid substrates, but is supplemented by some extra active terms (here the active stress σ^a and the active molecular field h^a) that break the TRS. This type of diffuse-interface model has been solved numerically as a minimal model for cell motility.^49,104

5 Applications 3: dynamics of thin active droplets on solid substrates

In this section, we consider a thin droplet of active polar fluids moving on a solid substrate in a simple two dimensions (xz plane), as schematically shown in Fig. 3b. The detailed dynamics of such a droplet can be described by the full dynamic equation system that is derived in the previous Section 4. However, here we present a thin film description of the active droplet on solid substrates where the lubrication approximation applies.^{101,102,105–107} In this case, the shape of the droplet is described by its thickness h(x,t) in the vertical z direction as a function of horizontal position x at time t, as shown in Fig. 3b. The hydrodynamic velocity field and the average orientation of the active polar filaments inside the two-dimensional drop are represented by v(x,z,t) = (u,w), and the polarization vector p(x,z,t), respectively. Similar problems have been studied by several groups^108–110 for thin films of passive liquid crystals, by Kitavtsev et al.¹¹¹ for thin films of active liquid crystals, by Sankararaman & Ramaswamy¹⁰⁵ for thin films of active polar fluids, and by Joanny & Ramaswamy¹⁰² for thin drops of active polar fluids. Here we use OVP to provide an alternative derivation of thin film equations for active polar droplets and to find approximate scaling laws for the spreading or dewetting of the droplet on solid substrate.

5.1 Thin film equations for active polar fluids

We firstly use OVP to derive the thin film equations for a thin active polar droplet moving on a solid substrate. To be specific and simple, we make several assumptions as follows.

(i) The active fluid is incompressible, satisfying the incompressibility condition, ∇·v = ∂_xu + ∂_zw = 0, from which we get the local mass conservation equation for the evolution of film height as

Here and in the following, we use ∂α to denote the partial derivatives with respect to the variable α such as time t, coordinates x and z.

(ii) The lubrication approximation^{66,77,112,113} is applied to the thin-film dynamics of the active polar droplet on the solid substrate. In the long-wave limit, the characteristic film thickness h₀ is much smaller than the length scale R₀ for variations in the x direction, i.e., h₀/R₀ ≪ 1. It follows that the film thickness varies slowly in space with |∂_xh| ≪ 1. Given h₀/R₀ ≪ 1 and ∇·v = 0, we obtain that the flow velocity v is approximately along the x direction with w ≪ u.

(iii) The equilibrium contact angle of the droplet θ_e is very small such that Young's equation is approximated as

where γ is fluid–gas surface tension, γ_SG solid–gas surface tension, and γ_SL solid–fluid surface tension.

(iv) We only consider droplet dynamics with left-right symmetry and the droplet shape is mainly determined by its interfacial energy and the effects of nematic elastic energy can be neglected. This arises when the characteristic thickness of the droplet is much larger than h_K ∼ K_p/γ with K_p being the elastic constant defined in eqn (32b). Then the total energy functional of the droplet is given by

where R is half of the contact length of the droplet with the solid substrate. The rate of change of the total energy is easily obtained aswhere we have used the identity ∂_th + Ṙ∂_xh = 0 (obtained from h(x = R(t),t) = 0) at the contact line x = R.

(v) The active filaments inside the droplet lack a head–tail polarity, that is, p and −p are equivalent, but they can show average nematic alignment. Furthermore, in the case of thin active droplets, the z dependence of p = (p_x,p_z) is determined by the equilibrium equations,

∂_z²p_x = ∂_z²p_z = 0, (61)

which are obtained from the minimization of nematic elastic energy, , similarly as defined in eqn (32b) under the approximation of one elastic constant, K_p, and fixed magnitude of p.

(vi) We assume the planar anchoring conditions at any bounding surfaces with which the active filaments are in contact, that is, the polarization vector p is parallel to the tangent direction of all the bounding surfaces. Here we then have: p = x̂ at z = 0 and p = [small tau, Greek, circumflex] ≈ (1,∂_xh) at the free surface z = h, in which x̂ is the unit vector along the x-direction and [small tau, Greek, circumflex] is the unit tangent vector of the free surface of the droplet. This anchoring boundary condition is mainly motivated by the stress-fiber structure in adherent cells² and by the experimental observations¹¹⁴ on thin films of amoeboid cells, in which the cells lie in the plane of the glass slide on which they spread and form nematic liquid-crystal structures. The planar anchoring conditions at the free surface have been employed in many previous works.^{98,102,105,106} In contrast, in our diffuse-interface model of active polar droplets in Section 4, we have assumed planar anchoring condition at the solid surface but perpendicular anchoring condition at the free interface, which mimics the orientation of actin filaments in the lamellipodium of migrating cells. Such anchoring boundary conditions have also been used in many previous works.^{43,49,101,104}

Then based on the assumptions in (v) and the planar anchoring conditions in (vi), we can solve the polarization vector p from eqn (61) independently of the flow velocity v for a given drop profile h(x,t) and obtain^102,105

p_x ≈ 1, p_z ≈ (z/h)∂_xh. (62)

(vii) We consider only the dynamic limit of active polar fluids discussed in Section 4.2, at which the effects of polarization and flow on the ATP hydrolysis are negligibly small, and the rate r ≈ Λ⁻¹Δμ is a constant parameter. In this case, the active stress σ^a = −ζpp in eqn (53a) is only a function of the local polarization vector p and it breaks the TRS, driving the system out of equilibrium locally. The work power done by σ^a on the thin droplet is approximated to the leading order as

(viii) Using the lubrication approximation, the dissipation functional is given to the leading order by^{66,77,112,113}

Here we have neglected the dissipation from the fluid slip at the solid surface (away from the contact line). The first term in Φ represents the usual hydrodynamic viscous dissipation under the lubrication approximation with η being the shear viscosity of the fluid. The second term represents the extra energy dissipation associated with the dynamics near the contact line region,^77,112 which can be very complex for the suspension droplets of active filaments.¹⁰¹ The phenomenological parameter ζ_cl is the frictional coefficient at the contact line, which is infinitely large for a pinned contact line, and is zero for a freely moving contact line.

From the above discussions, we then obtain the Rayleighian as

where p is the pressure (a Lagrange multiplier) that imposes the incompressibility constraint. Minimizing [scr R, script letter R] with respect to the rates, u, w, ∂_th, and Ṙ, give a closed dynamic equation system for the thin droplets moving on solid substrates as follows. In the bulk, the dynamic equations are

∂_zp = 0, (66b)

which are supplemented with the no-slip condition u = 0 at the solid surface z = 0, the boundary conditions at the free surface z = h:

p(x,h,t) = p₀ − γ∂_x²h, (67a)

η∂_zu|_z=h = 0, (67b)

and the boundary conditions at the contact line x = R and z = 0:Here p₀ is the pressure in the surrounding gas; we have also used the impermeability condition w = 0 at the solid surface z = 0, and the kinematic boundary condition w = ∂_th + u∂_xh at the free surface z = h.

As in the classical problems of thin film fluids,^66,107,115 the solution of the above closed equation system gives a parabolic profile for u(x,z,t) in the form of

from which we find the thickness-averaged velocity asSubstituting eqn (70) into the mass conservation eqn (57), one obtains the thin film equation for active droplets:which is supplemented with the boundary condition (68) at the contact line x = R(t). The static solutions of this thin film equation yield the steady-state shape of the active droplet and the dynamic scaling properties of the solutions lead to the spreading or dewetting laws for the active droplet.^32,33

We would like to comment and compare our model for the thin active droplets on solid substrates with other models in the literature^101,102,106 as follows.

(i) In comparison to the thin-film model by Joanny & Ramaswamy,¹⁰² we have neglected the effect of nematic elastic energy in determining the droplet shape. As a result, our thin film equation is a limiting case of their model when the droplet thickness is much larger than h_K ∼ K_p/γ as discussed above near eqn (59). However, if h_K is not very small, our boundary condition (68) in the vicinity of the contact line with h_K > h ∼ 0 may be problematic and elastic contributions have to be included.

(ii) In the thin-film model by Loisy et al.,¹⁰⁶ the flows inside the active droplets are induced by the winding of the polarization field. This winding introduces a dramatic change in the orientation of the polarization p along the thickness z-direction: , in which θ is the angle of p relative to the x-axis and ω is an integer winding number that counts the number of quarter turns of p across the drop height. In comparison, in the present work, we have not considered such internal polarization winding and the orientation of p varies along z as θ ∼ p_z/p_x ≈ (z/h)∂_xh (see eqn (62)). Such difference in the variation of p orientation leads to the difference in the final form of thin film equation between the present work and that by Loisy et al.¹⁰⁶

(iii) In the recent work by Trinschek et al.,¹⁰¹ the authors have proposed a more complete model for active polar droplets, which is similar to our model presented in Section 4 but have introduced one additional active contribution from the treadmilling or self-propulsion of active units in the direction of their polarization. Using their more complete free energy, we can still apply our variational approach and the lubrication approximation to study the more complicated thin film dynamics by following similar methods that we have done for thin films of binary mixtures before.⁶⁶

5.2 Spreading laws for thin active fluid droplets

Now we consider the spreading dynamics of thin active droplets on solid substrates. We do not try to solve the thin film eqn (71) directly, but use OVP as an approximation method to solve the scaling laws for the droplet spreading.

We assume that the height profile of the droplet h(x,t) is given by a parabolic function

and the velocity u(x,z,t) inside the droplet takes the following parabolic profile along the z-direction asThe reason of choosing these function forms is apparent if we compare eqn (72b) with the velocity profile in eqn (69), and remember that the equilibrium shape of a thin droplet takes the form of eqn (72a) with θ(t) = θ_e being the equilibrium contact angle and R(t) = R_e being (half of) the contact length between the droplet and the solid substrate at equilibrium. Note that θ(t) defined from the parabolic height profile in eqn (72a) is the apparent contact angle interpolated away from the microscopic contact line, which is different from the contact angle θ = |∂_xh| defined in the previous subsection locally in the close vicinity of the contact line.

To achieve an approximate description of the droplet dynamics, the time-dependent parameters θ(t), R(t), and u₀(x,t) must be determined by OVP. However, note that these parameters are not independent. Firstly, from the conservation of the droplet area (or mass), A₀, that is, , we have

Secondly, from the mass conservation eqn (57) for droplet height, we obtainSubstituting the approximate profile of h(x,t) and u(x,z,t) in eqn (72) into eqn (74) we obtainTherefore, the dynamics of the droplet can be described by one time-dependent parameter and here we take R(t).

Substituting the droplet profile h(x,t) in eqn (72a) into eqn (59), we obtain the total free energy

from which we find the rate of change of the total energy asSimilarly, we find the energy dissipation function Φ asand the work power done by active stress, σ^a, aswhich are obtained by substituting the approximate profile of h(x,t) and u(x,z,t) in eqn (72) into the dissipation function in eqn (64) and the active work power function in eqn (63), respectively. Here C ≡ 2[ln(2R/ε) − 2] and ε is the molecular cutoff length that is introduced to remove the divergence in the energy dissipation near the contact line.

Then from eqn (77)–(79), we obtain the Rayleighian . Minimizing [scr R, script letter R] with respect to Ṙ gives the following evolution equation

Here sign(ζ) is the sign function of ζ with sign(ζ) = +1 for extensile active units and sign(ζ) = −1 for contractile active units as shown in Fig. 2. The contact angle θ is a function of R as given by eqn (73). The dimensionless parameter k_cl ≡ ζ_cl/ζ_hydro is a material parameter determined by the droplet and the substrate, which can be treated as a constant and characterizes the importance of the additional friction ζ_cl near the contact line relative to the normal hydrodynamic friction ζ_hydro = 3Cη/2θ in the bulk fluids of the droplet. Two time scales are introduced and defined asThe time τ_act represents the characteristic time for the droplet to reach steady-state motion driven by active stresses, and τ_rel represents the relaxation time needed for the droplet to reach the equilibrium contact angle θ_e. The dimensionless parameter k_act defined bycharacterizes the strength of activity in the active fluids. If k_act is small, the equilibration of the droplet shape and contact angle is very fast and the droplet spreading is mainly driven by surface energy. On the other hand, if k_act is large, the activity is strong and the activity-driven droplet motion is much faster than the energy-driven equilibration of the droplet, that is, the droplet spreading is mainly driven by active stress and a large capillary-number flow will be induced.

Particularly, for k_act ≪ 1, the first term on the right-hand side of eqn (80) can be ignored, and the evolution equation for R becomes

which gives (for θ_e → 0) the classical Tanner's spreading law for two-dimensional droplets:^32–34On the other hand, if k_act ≫ 1, the second term on the right-hand side of eqn (80) can be ignored, and the equation becomesfrom which we obtain the spreading law (or dewetting law from very small initial contact angle) predicted by Joanny & Ramaswamy¹⁰² for two-dimensional droplets of extensile active units with sign(ζ) = +1 (or of contractile active units with sign(ζ) = −1):

Note that as mentioned by Joanny & Ramaswamy,¹⁰² the effects of activity on the droplet spreading enter at the same order in gradients as those of gravity, but with a different dependence on the film height. Furthermore, similar dynamic equation as eqn (80) for thin droplets on solid substrates has been obtained in a very different scenario where evaporation occurs at the free surface of the droplet.⁷⁷ The effects of activity on the droplet spreading enter at the same order as those of evaporation, but with a different dependence on the droplet radius or contact length.

In addition, the formulation and calculations presented here can be readily extended to the thin-film dynamics of three-dimensional droplets on solid substrates, particularly for the spreading dynamics of a droplet with cylindrical symmetry. Furthermore, the effects of nematic energy on the spreading dynamics can also be considered by including nematic elastic energy,¹⁰² which takes the simple form of .

6 Conclusions

Onsagers variational principle (OVP) has recently become an indispensable and powerful tool in the study of the nonlinear and nonequilibrium phenomena of many inert soft matter systems, such as liquid droplets, colloid suspensions, nematic liquid crystals, polymer gels, and surfactants, etc. In this work, we present a simple extension of OVP for the dynamic modeling of active soft matter such as suspensions of bacteria and aggregates of animal cells. In this extended OVP, the active forces generated locally by individual active units are treated as non-conservative forces that cannot be derived from any free energy and dissipation functions. We then apply this extended form of OVP to three representative active matter problems, which are motivated by the biology of bacteria and animal cells. We show that OVP can not only help to formulate thermodynamically-consistent models, but can also be used to find approximate solutions for the emergent structures and complex dynamics of active soft matter.

The first application of OVP presented in Section 3 is about the directional motion of individual active units: a molecular motor walking on a stiff biofilament and a toy two-sphere microswimmer moving in a viscous fluid. In the motor/filament system, we consider the mechanochemical cross-coupling which indicates that the system is a reversible machine: it can not only convert chemical energy into mechanical work, but can also convert mechanical work into chemical energy. In the toy microswimmer, we show how directional self-propulsion can be generated by cyclic body-shape oscillations together with front-back asymmetry in hydrodynamic friction. It is shown that mechanochemical cross-coupling in biological systems can be considered in Onsager's framework of non-equilibrium thermodynamics. Activity and the broken time reversal symmetry in active matter are basically resulted from the persistent consumption and conversion of chemical energy, released during spontaneous ATP hydrolysis, into motion or mechanical work.

The second application presented in Section 4 is about the two-phase hydrodynamics for a droplet of active polar fluids, which is composed of suspending contractile or extensile active units such as bacteria, actomyosin units, and animal cells. We use OVP to formulate a diffuse-interface model for an active polar droplet moving on a solid substrate. This hydrodynamic model is thermodynamically consistent in both hydrodynamic equations in the bulk fluid and matching boundary conditions at the solid surface.

The third application presented in Section 5 is about the motion of a thin active polar droplet on a solid substrate in two dimensions. Using the lubrication approximation, we firstly apply OVP to derive the classical thin film equation that has been obtained previously. We then use OVP as an approximation tool to find two scaling laws for the spreading (or dewetting) of the thin active droplet in the respective limits of negligible activity and strong activity. It is interesting to note that the reduced equation obtained for the spreading (or dewetting) dynamics of thin active droplets takes a similar form as that has been obtained previously for the dewetting dynamics of an evaporating droplet on solid substrates.

Below we make a few general remarks and outlook.

(i) Near-equilibrium assumption of OVP. OVP is proposed in Onsager's linear-response framework of non-equilibrium thermodynamics, which is based on the near-equilibrium assumption. However, biological systems are usually far away from equilibrium. Therefore, the validity and the range of the OVP applications should be and can only be justified by solving real biological problems and comparing with quantitative experiments.^1,4

(ii) Relationships between OVP/OMVP and other approaches in nonequilibrium thermodynamics. Following the pioneering works of Onsager, there have been various approaches developed for nonequilibrium thermodynamics. In particular, there have been various variational principles formulated for the study of irreversible processes.¹¹⁶ More recently, the general equation for non-equilibrium reversible-irreversible coupling (GENERIC) formalism has been proposed as an extension of Hamiltonian's formalism of classical mechanics to nonequilibrium thermodynamic systems with both reversible and irreversible dynamics.¹¹⁷ However, a general discussion on the relationships among the various approaches is beyond the scope of this work.

(iii) Applications of Onsager–Machlup variational principle (OMVP). In the end of Section 2.1, we have mentioned that Onsager and Machlup⁶⁰ introduced OMVP in their study of the statistical fluctuations of kinetic paths in the framework of Langevin equation. They have shown that the most probable kinetic path over a certain long-time period is determined by the minimization of a time integral, i.e., the Onsager–Machlup integral in eqn (12). Recently, Doi et al.⁶⁸ proposed that OMVP can be used to approximate the long-time dynamics of nonequilibrium systems. However, in this work, we have not given applications of using OMVP to find approximate solutions for long-time behaviors such as steady states. We mention two potential applications of OMVP as follows: the steady-state for the wave propagation and sustained oscillations observed in migrating cells;^29–31 the steady-state (spontaneous) retraction dynamics of an injured axon¹¹⁸ or a laser-cutting stress-fiber bundle in adherent cells.¹¹⁹

(iv) Applications of OVP and OMVP to more specific biological problems. The applications considered in this work are mostly toy models or simplified models of mostly theoretical interests. We are now trying to apply the extended form of OVP to more specific biological problems such as cell spreading, cell curvotaxis, wound closure, tissue folding, and so on. However, in these real systems, we usually need to involve many more complex active processes¹⁶ in addition to active forces and cyclic body-shape oscillations, such as tensional homeostasis,² cell division and apoptosis,¹³ topological cell rearrangements,¹⁶ memory effects.^1,17

In summary, the variational method proposed in this work about incorporating biochemical activity into OVP will help to construct thermodynamically-consistent models and to find approximate dynamic solutions in active soft matter. Particularly, this will help to deepen our understanding of the emergent structure and dynamic behaviors of real in vivo biological systems such as bacteria suspensions, individual animal cells and cell aggregates (or tissues).^1,2,13

Conflicts of interest

There are no conflicts to declare.

A Some additional notes on the choice of thermodynamic fluxes and forces

In Onsager's theory of non-equilibrium thermodynamics, the choice of thermodynamic fluxes and forces is not unique.^50,51 For example, for active polar fluids, two different sets of fluxes and forces can be chosen as follows:

(i) Fluxes are chosen to be r, P, [small epsi, Greek, dot above] (or the flow velocity v) and the corresponding forces are Δμ, h, σ^v, respectively, as taken in Section 4 of this work. Here all the three fluxes have the same time parity.

(ii) Fluxes are chosen to be r, P, σ^v (the momentum flux) and the corresponding forces are Δμ, h, [small epsi, Greek, dot above], respectively, as taken by Marchetti et al.¹ Here the time parity of the flux σ^v is different from the other two fluxes r and P. Furthermore, note that in this case, the new pair of flux σ^v and force [small epsi, Greek, dot above] is a swap of the pair of flux [small epsi, Greek, dot above] and force σ^v in the first choice (i).

The symmetry of Onsager matrix coupling fluxes and forces depends on the time parity (i.e., the time-reversal signature) of the fluxes:^1,50 the Onsager coupling matrix is symmetric for fluxes of the same time parity and is antisymmetric for fluxes of opposite time parity. Therefore, for the first choice (i) with fluxes of the same time parity, the Onsager coupling matrix is symmetric. In comparison, for the second choice (ii), the cross-coupling coefficients between the flux σ^v and the flux r or between the flux σ^v and the flux P are both antisymmetric as shown in Marchetti et al.¹

Given the non-unique choice of flux–force pairs, we will make some general discussions in this Appendix about the consequences of different choices of thermodynamic fluxes and forces.

A.1 Changes in the symmetry of Onsager coupling matrix by swapping some of fluxes and forces

We discuss the changes in the symmetry of Onsager coupling matrix when we swap some of fluxes with forces, as in the above two choices of force–flux pairs in active polar fluids. To this end, let's consider two pairs of fluxes and forces: [small alpha, Greek, dot above]₁ and X₁, [small alpha, Greek, dot above]₂ and X₂. Suppose that the two fluxes [small alpha, Greek, dot above]₁ and [small alpha, Greek, dot above]₂ have the same time parity, and hence the two forces X₁ and X₂ also have the same time parity but opposite to that of the corresponding fluxes. In this case, the linear force–flux relations are then given by:in which the Onsager coupling matrix L_ij is symmetric and positive definiteness, that is, L₁₂ = L₂₁, L₁₁, L₂₂ > 0, and L₁₁L₂₂ − L₁₂² > 0.

If alternatively, we swap X₂ with [small alpha, Greek, dot above]₂, that is, we choose the fluxes to be [small alpha, Greek, dot above]₁ and X₂, and the corresponding forces are X₁ and [small alpha, Greek, dot above]₂. Then the time parities of the two fluxes [small alpha, Greek, dot above]₁ and X₂ are now different. In this case, the linear force–flux relations become

in which L̃₁₁ = (L₁₁L₂₂ − L₁₂²)/L₂₂, L̃₁₂ = L₁₂/L₂₂, L̃₂₁ = −L₂₁/L₂₂, L̃₂₂ = 1/L₂₂. Note that now L̃₂₁ = −L̃₁₂, that is, the Onsager coupling matrix L̃_ij is antisymmetric. From this simple example, we show that although there exists some flexibility in choosing fluxes and forces, we can safely use the time parity of the chosen fluxes to determine the symmetry of the Onsager coupling matrix in the linear force–flux relations.

A.2 Onsager's variational principle for thermodynamic fluxes of different time parities

Suppose that there are two sets of slow state variables, α and β and they have opposite time parity. Then the free energy is [scr F, script letter F](α,β) and the rate of change of the free energy is given byin which [small alpha, Greek, dot above] and are two fluxes, and

X_i ≡ − ∂[scr F, script letter F]/∂α_i, Y_i ≡ − ∂[scr F, script letter F]/∂β_i (A4)

are the corresponding conjugate forces. Since α and β have opposite time parities, then the time parities of [small alpha, Greek, dot above] and , and that of X and Y are both opposite, respectively. In this case, we can decompose the two forces into reactive and dissipative parts as

X = X^d + X^r, Y = Y^d + Y^r. (A5)

Here X^d and Y^d arise from dissipative couplings that can be derived from dissipation function, while X^r and Y^r derive from reactive couplings:with the friction matrix L^αβ_ij = −L^βα_ji being antisymmetric. The reactive couplings do not contribute to dissipation and therefore can not be derived from dissipation functions. Substituting eqn (A5) into the rate of change of the free energy in eqn (A3), using eqn (A6) and L^αβ_ij = −L^βα_ji, we obtain

The dissipation function is given by

with the friction matrix L^αα_ij = L^αα_ji and L^ββ_ij = L^ββ_ji being symmetric. Using eqn (A5)–(A8) and minimizing the Rayleighian , we obtainor equivalently by taking the inverse:where the mobility matrix R is the inverse of the friction matrix L. Note that only the dissipative part of the thermodynamic forces can be considered and derived from the minimization of the Rayleighian, which is simply because the reactive forces don't contribute to the dissipation function and the Rayleighian.

Acknowledgements

We thank Masao Doi at Beihang University for many fruitful discussions, particularly during his visit in Xu's group in 2019. Great thanks should also be given to Len Pismen and Yariv Kafri at Technion for their valuable comments and suggestions. T. Q. acknowledges support of the Hong Kong RGC CRF (No. C1018-17G). X. X. is supported in part by National Natural Science Foundation of China (NSFC, No. 12004082), by Provincial Science Foundation of Guangdong (No. 2019A1515110809), by Guangdong Province Universities and Colleges Pearl River Scholar Funded Scheme (2019), by 2020 Li Ka Shing Foundation Cross-Disciplinary Research Grant (No. 2020LKSFG08A), by Guangdong Basic and Applied Basic Research Foundation (No. 2020B1515310005), and by Featured Innovative Projects (No. 2018KTSCX282), Youth Talent Innovative Platforms (No. 2018KQNCX318) in Universities in Guangdong Province and Science and Technology Innovation Program 2020 (No. ST002).

Notes and references

1. M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189.
2. U. S. Schwarz and S. A. Safran, Rev. Mod. Phys., 2013, 85, 1327–1381.
3. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
4. J. Prost, F. Jülicher and J.-F. Joanny, Nat. Phys., 2015, 11, 111–117.
5. A. M. Menzel, Phys. Rep., 2015, 554, 1–45.
6. G. De Magistris and D. Marenduzzo, Phys. A, 2015, 418, 65–77.
7. D. Marenduzzo, Eur. Phys. J.: Spec. Top., 2016, 225, 2065–2077.
8. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
9. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601.
10. J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
11. J. Gachelin, A. Rousselet, A. Lindner and E. Clement, New J. Phys., 2014, 16, 025003.
12. H.-P. Zhang, A. Beer, E.-L. Florin and H. L. Swinney, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 13626–13630.
13. J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost and F. Jülicher, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 20863–20868.
14. M. H. Köpf and L. M. Pismen, Soft Matter, 2013, 9, 3727–3734.
15. S. He, Y. Green, N. Saeidi, X. Li, J. J. Fredberg, B. Ji and L. M. Pismen, J. Mech. Phys. Solids, 2020, 137, 103860.
16. M. Popović, A. Nandi, M. Merkel, R. Etournay, S. Eaton, F. Jülicher and G. Salbreux, New J. Phys., 2017, 19, 033006.
17. S. Banerjee and M. C. Marchetti, in Cell Migrations: Causes and Functions, ed. C. A. M. La Porta and S. Zapperi, Springer, Cham, 2019, pp. 45–66.
18. J. Toner and Y. Tu, Phys. Rev. Lett., 1995, 75, 4326–4329.
19. C. K. Hemelrijk and H. Hildenbrandt, Interface Focus, 2012, 2, 726–737.
20. D. Helbing, Rev. Mod. Phys., 2001, 73, 1067.
21. C. Castellano, S. Fortunato and V. Loreto, Rev. Mod. Phys., 2009, 81, 591.
22. I. S. Aranson, D. Volfson and L. S. Tsimring, Phys. Rev. E, 2007, 75, 051301.
23. C. Scholz, M. Engel and T. Pöschel, Nat. Commun., 2018, 9, 1–8.
24. J. Zhang, E. Luijten, B. A. Grzybowski and S. Granick, Chem. Soc. Rev., 2017, 46, 5551–5569.
25. Z. Lin, C. Gao, M. Chen, X. Lin and Q. He, Curr. Opin. Colloid Interface Sci., 2018, 35, 51–58.
26. Z. Ma, M. Yang and R. Ni, Adv. Theory Simul., 2020, 3, 2000021.
27. A. Grosberg and J.-F. Joanny, Phys. Rev. E, 2015, 92, 032118.
28. M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys., 1993, 65, 851–1112.
29. G. Gerisch, T. Bretschneider, A. Müller-Taubenberger, E. Simmeth, M. Ecke, S. Diez and K. Anderson, Biophys. J., 2004, 87, 3493–3503.
30. O. D. Weiner, W. A. Marganski, L. F. Wu, S. J. Altschuler and M. W. Kirschner, PLoS Biol., 2007, 5, e221.
31. N. Inagaki and H. Katsuno, Trends Cell Biol., 2017, 27, 515–526.
32. L. Leger and J.-F. Joanny, Rep. Prog. Phys., 1992, 55, 431.
33. P.-G. De Gennes, Rev. Mod. Phys., 1985, 57, 827.
34. D. Bonn, J. Eggers, J. Indekeu, J. Meunier and E. Rolley, Rev. Mod. Phys., 2009, 81, 739.
35. H. Tanaka, Faraday Discuss., 2012, 158, 371–406.
36. R. G. Winkler, D. A. Fedosov and G. Gompper, Curr. Opin. Colloid Interface Sci., 2014, 19, 594–610.
37. P.-G. De Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, 1993.
38. S. B. Ross-Murphy, Polym. Gels Networks, 1994, 2, 229–237.
39. M. Doi, J. Phys. Soc. Jpn., 2009, 78, 052001.
40. M. Cates and S. Candau, J. Phys.: Condens. Matter, 1990, 2, 6869–6892.
41. S. Safran, P. Pincus, D. Andelman and F. MacKintosh, Phys. Rev. A, 1991, 43, 1071–1078.
42. A. Baskaran and M. C. Marchetti, J. Stat. Mech.: Theory Exp., 2010, 2010, P04019.
43. M. E. Cates and E. Tjhung, J. Fluid Mech., 2018, 836, P1.
44. R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101.
45. J. Prost and R. Bruinsma, Europhys. Lett., 1996, 33, 321–326.
46. E. Woillez, Y. Zhao, Y. Kafri, V. Lecomte and J. Tailleur, Phys. Rev. Lett., 2019, 122, 258001.
47. T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen and O. Shochet, Phys. Rev. Lett., 1995, 75, 1226.
48. R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo and M. E. Cates, Nat. Commun., 2014, 5, 1–9.
49. E. Tjhung, A. Tiribocchi, D. Marenduzzo and M. Cates, Nat. Commun., 2015, 6, 1–9.
50. S. R. de Groot and P. Mazur, Non-equilibrium Thermodynamics, Dover, New York, 1984.
51. I. Gyarmati, Non-equilibrium Thermodynamics, Springer, Berlin, 1970.
52. H. Li, X.-Q. Shi, M. Huang, X. Chen, M. Xiao, C. Liu, H. Chaté and H. Zhang, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 777–785.
53. M. Doi, Soft Matter Physics, Oxford University Press, 2013.
54. L. Onsager, Phys. Rev., 1931, 37, 405–426.
55. L. Onsager, Phys. Rev., 1931, 38, 2265–2279.
56. J. W. Strutt, Proc. London Math. Soc., 1871, 1, 357–368.
57. L. Onsager, Ann. N. Y. Acad. Sci., 1945, 46, 241–265.
58. T. Qian, X.-P. Wang and P. Sheng, J. Fluid Mech., 2006, 564, 333–360.
59. M. Doi and S. Edwards, The Theory of Polymer Dynamics, Oxford University Press, 1986.
60. L. Onsager and S. Machlup, Phys. Rev., 1953, 91, 1505–1512.
61. S. Ono, Adv. Chem. Phys., 1961, 3, 267–321.
62. M. Doi, J. Condens. Matter Phys., 2011, 23, 284118.
63. M. Doi, Prog. Polym. Sci., 2021, 112, 101339.
64. Y. Wang, C. Liu, P. Liu and B. Eisenberg, 2020, arXiv:2001.10149.
65. Q. Wang, in Frontiers and Progress of Current Soft Matter Research, ed. X.-Y. Liu, Springer, Singapore, 2020, pp. 101–132.
66. X. Xu, U. Thiele and T. Qian, J. Phys.: Condens. Matter, 2015, 27, 085005.
67. Y. Di, X. Xu, J. Zhou and M. Doi, Chin. Phys. B, 2018, 27, 024501.
68. M. Doi, J. Zhou, Y. Di and X. Xu, Phys. Rev. E, 2019, 99, 063303.
69. X. Xu and T. Qian, Procedia IUTAM, 2017, 20, 144–151.
70. P. Sheng, J. Zhang and C. Liu, Prog. Theor. Phys. Suppl., 2008, 175, 131–143.
71. T. V. Sachin Krishnan, R. Okamoto and S. Komura, Phys. Rev. E, 2016, 94, 062414.
72. Y. Oya and T. Kawakatsu, J. Chem. Phys., 2018, 148, 114905.
73. M. Arroyo, N. Walani, A. Torres-Sánchez and D. Kaurin, The Role of Mechanics in the Study of Lipid Bilayers, Springer, Cham, 2018, pp. 287–332.
74. A. Fang, T. Qian and P. Sheng, Phys. Rev. E, 2008, 78, 061703.
75. X. Xu, X. Man, M. Doi, Z.-C. Ou-Yang and D. Andelman, Macromolecules, 2019, 52, 9321–9333.
76. M. Doi, Chin. Phys. B, 2015, 24, 20505.
77. X. Man and M. Doi, Phys. Rev. Lett., 2016, 116, 066101.
78. X. Man and M. Doi, Phys. Rev. Lett., 2017, 119, 044502.
79. J. Zhou and M. Doi, Phys. Rev. Fluids, 2018, 3, 084004.
80. J. Sui, M. Doi and Y. Ding, Soft Matter, 2018, 14, 8956–8961.
81. P. Khunpetch, X. Man, T. Kawakatsu and M. Doi, J. Chem. Phys., 2018, 148, 134901.
82. F. J. Vernerey and U. Akalp, Phys. Rev. E, 2016, 94, 012403.
83. Y.-H. Zhang, M. Deserno and Z.-C. Tu, Phys. Rev. E, 2020, 102, 012607.
84. X. Yang, J. Li, M. G. Forest and Q. Wang, Entropy, 2016, 18, 202.
85. A. Tiribocchi, R. Wittkowski, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2015, 115, 188302.
86. S. Gu, T. Qian, H. Zhang and X. Zhou, Chaos, 2020, 30, 053133.
87. Q. Du, T. Li, X. Li and W. Ren, Sci. China Math., 2020, 1–42.
88. H. Touchette, Phys. Rep., 2009, 478, 1–69.
89. J. N. Reddy, Energy Principles and Variational Methods in Applied Mechanics, John Wiley & Sons, 2017.
90. B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts and P. Walter, Molecular Biology of the Cell, Garland Science, London, 5th edn, 2007.
91. O. Kedem and S. R. Caplan, Trans. Faraday Soc., 1965, 61, 1897–1911.
92. Y.-D. Chen and T. L. Hill, Proc. Natl. Acad. Sci. U. S. A., 1974, 71, 1982–1986.
93. A. Patel, L. Malinovska, S. Saha, J. Wang, S. Alberti, Y. Krishnan and A. A. Hyman, Science, 2017, 356, 753–756.
94. R. Phillips, J. Kondev, J. Theriot and H. Garcia, Physical Biology of the Cell, Garland Science, 2012.
95. F. Jülicher, A. Ajdari and J. Prost, Rev. Mod. Phys., 1997, 69, 1269–1281.
96. H. C. Astley, C. Gong, J. Dai, M. Travers, M. M. Serrano, P. A. Vela, H. Choset, J. R. Mendelson, D. L. Hu and D. I. Goldman, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 6200–6205.
97. D. L. Bodor, W. Pönisch, R. G. Endres and E. K. Paluch, Dev. Cell, 2020, 52, 550–562.
98. A. Loisy, J. Eggers and T. B. Liverpool, Soft Matter, 2020, 16, 3106–3124.
99. K. Zimmermann, I. Zeidis and C. Behn, Mechanics of Terrestrial Locomotion: with a Focus on Non-Pedal Motion Systems, Springer, 2009.
100. A. Lurie, Analytical Mechanics, Springer, 2002.
101. S. Trinschek, F. Stegemerten, K. John and U. Thiele, Phys. Rev. E, 2020, 101, 062802.
102. J.-F. Joanny and S. Ramaswamy, J. Fluid Mech., 2012, 705, 46.
103. A. C. Callan-Jones and F. Jülicher, New J. Phys., 2011, 13, 093027.
104. E. Tjhung, D. Marenduzzo and M. E. Cates, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 12381–12386.
105. S. Sankararaman and S. Ramaswamy, Phys. Rev. Lett., 2009, 102, 118107.
106. A. Loisy, J. Eggers and T. B. Liverpool, Phys. Rev. Lett., 2019, 123, 248006.
107. W. Ren, D. Hu and W. E, Phys. Fluids, 2010, 22, 102103.
108. M. Ben Amar and L. Cummings, Phys. Fluids, 2001, 13, 1160–1166.
109. S. Mechkov, A.-M. Cazabat and G. Oshanin, J. Condens. Matter Phys., 2009, 21, 464131.
110. T.-S. Lin, L. J. Cummings, A. J. Archer, L. Kondic and U. Thiele, Phys. Fluids, 2013, 25, 082102.
111. G. Kitavtsev, A. Münch and B. Wagner, Proc. R. Soc. A, 2018, 474, 20170828.
112. H. Wang, D. Yan and T. Qian, J. Phys.: Condens. Matter, 2018, 30, 435001.
113. R. V. Craster and O. K. Matar, Rev. Mod. Phys., 2009, 81, 1131.
114. R. Kemkemer, V. Teichgräber, S. Schrank-Kaufmann, D. Kaufmann and H. Gruler, Eur. Phys. J. E, 2000, 3, 101–110.
115. A. Oron, S. H. Davis and S. G. Bankoff, Rev. Mod. Phys., 1997, 69, 931–980.
116. M. Ichiyanagi, Phys. Rep., 1994, 243, 125–182.
117. H. C. Öttinger, Beyond equilibrium thermodynamics, John Wiley & Sons, 2005.
118. X. Shao, R. You, T. H. Hui, C. Fang, Z. Gong, Z. Yan, R. C. C. Chang, V. B. Shenoy and Y. Lin, Biophys. J., 2019, 117, 193–202.
119. A. Besser, J. Colombelli, E. H. Stelzer and U. S. Schwarz, Phys. Rev. E, 2011, 83, 051902.

---

Footnotes

† It is interesting to mention that in active matter, the presence of active forces and the breakdown of TRS can significantly change the statistical behaviors of the system such as the distribution of α and barrier crossing kinetics. The dynamics of downhill and uphill processes are very different even for passive systems, and hence must be treated using their respective variational approaches. In this work, we are interested in the active dynamics that would reduce to downhill processes if the activity vanishes. However, the fluctuation effects can still be investigated in a variational framework.^46,86

‡ Private communications with M. Doi.

---