The role of fluctuations in determining cellular network thermodynamics.

The steady state distributions of phenotypic responses within an isogenic population of cells result from both deterministic and stochastic characteristics of biochemical networks. A biochemical network can be characterized by a multidimensional potential landscape based on the distribution of responses and a diffusion matrix of the correlated dynamic fluctuations between N-numbers of intracellular network variables. In this work, we develop a thermodynamic description of biological networks at the level of microscopic interactions between network variables. The Boltzmann H-function defines the rate of free energy dissipation of a network system and provides a framework for determining the heat associated with the nonequilibrium steady state and its network components. The magnitudes of the landscape gradients and the dynamic correlated fluctuations of network variables are experimentally accessible. We describe the use of Fokker-Planck dynamics to calculate housekeeping heat from the experimental data by a method that we refer to as Thermo-FP. The method provides insight into the composition of the network and the relative thermodynamic contributions from network components. We surmise that these thermodynamic quantities allow determination of the relative importance of network components to overall network control. We conjecture that there is an upper limit to the rate of dissipative heat produced by a biological system that is associated with system size or modularity, and we show that the dissipative heat has a lower bound.

Introduction

Measurements of individual cells within a population indicate that phenotypic differences between isogenic cells is common, even when they are in a homogeneous and stable environment. Despite differences between individuals, populations of cells can manifest apparently stable distributions of phenotypic expression. The measured phenotypic parameters can include concentration of protein products of gene expression, indicators associated with promoter activation, RNA transcripts, and complex cell traits such as morphology. The steady state distribution of phenotypes observed from single cell analysis is a probability density function and can be represented as a potential energy landscape [1-5]. Living cells are a clear example of a nonequilibrium system [6].

Studies that involve imaging of live cells reveal that despite the observed consistency of a steady state distribution of phenotypes across a population of cells, there can be significant dynamic variability in each cell and from one cell to another [3, 7–13]. A number of studies provide direct evidence that the distribution of phenotypes in steady state distributions is ergodic in that subsequent to a transient perturbation the population will eventually relax to the steady state distribution [3, 14–18] when culture conditions are kept constant. The recapitulation of the steady state distribution, even after single cell cloning, demonstrates that each cell or its progeny can explore all microstates on the landscape. The population response appears invariant, but the individual entities (cells) that comprise the population present a dynamic, random expression of phenotypes, which ultimately results in the heterogeneity of the population. These are characteristics of many physical systems that can be well-described by statistical mechanics. The variations in populations of cells are often attributed to stochastic fluctuations, or noise, due to small numbers of molecules associated with transcription and translation. However, fluctuations in small numbers of molecules is not the only source of this variation [9, 19–21]; correlation analysis has suggested that the main source of noise may be upstream regulatory components [14]. Ensembles of biochemicals are responsible for producing an observed phenotype [9, 21, 22], and these network components demonstrate coordinated concentration response functions. Importantly, dynamic fluctuations of network components, and correlations in fluctuations among multiple network components, are defining features of regulated networks [2, 13, 16, 21, 23–31]. An abbreviated biochemical network system is depicted in the schematic in Fig 1, based on the interactions between transcription factors involved in maintaining pluripotency [32], and provides an illustration of some of the general network features that our model addresses.

Fig 1

Schematic of an intracellular network.

This is based on concepts presented by Rizzino and Wuebben [32] on control of pluripotency, and consists of N = 3 variable components, or dimensions, of a network with activation and repression relationships as indicated. As described in the text, these relationships define the measurable covariances in an NxN diffusion matrix, which contribute to the thermodynamics of the network together with the gradients of a landscape (which are derived from the multi-dimensional steady state probability density). The nonequilibrium steady state in this open thermodynamic system is supported by an influx of free energy from outside the system, which is dissipated as heat.

Identification of network components and the nature of their interdependence remains a challenge. While some “omics” analysis methods can probe many variables simultaneously and in some cases at the level of individual cells, they provide only a snapshot in time. While these methods allow determination of the coincidence of molecular species, without dynamic data in individual cells that can reveal the rates and magnitudes of the interactions between network species, it is impossible to determine unambiguously the relationship between them or to determine a physical mechanistic basis for control of the network [11, 21, 33, 34]. A promising approach for observing single cell dynamics is imaging of live cells over time, which provides access to fluctuations in phenotypic expression in individual cells across a population [10–12, 16, 26]. Live cell imaging can also provide spatial and temporal information across scales and can provide ancillary information such as direct observation of the timing of cell division. We have previously used quantitative microscopy [3, 8] to track temporal responses of populations of individual cells, and Langevin/Fokker-Planck (L/FP) dynamics to analyze the population distributions as potential energy landscapes [3]. A cell line expressing green fluorescent protein (GFP) associated with the promoter for the gene for the extracellular matrix protein, tenascin C, was monitored by live cell fluorescence microscopy. The steady state distribution of expression levels of GFP in individual cells allowed construction of the potential landscape for the population. Changes in fluorescence intensity in individual cells were quantified at 15-minute intervals. By measuring real time trajectories of the activity of the gene promoter, we determined the diffusion coefficient for fluctuations in the promoter. We selected four subpopulations of cells by flow sorting, and by employing L/FP dynamics, accurately predicted the time-dependent relaxation of the subpopulations to the steady state distribution. The subpopulations were selected based on different levels of promoter activity and they relaxed to steady state with very different kinetics, but the diffusion coefficient provided excellent predictions of all relaxation kinetics with no adjustable parameters [3]. The study demonstrated that accurate experimental data on gene expression fluctuations can be collected by fluorescence microscopy, that data at the cellular level can provide additional details about control of gene expression that cannot be determined at the population level, and that FP dynamics is an appropriate approach to probing intracellular control mechanisms from experimental observations.

While such a one-dimensional landscape is a simple manifestation of a complex adaptive system, it provides little insight into the mechanisms responsible for establishing and maintaining the stability and dynamics of the steady state distribution because of the contributions from unseen variables. In the current theoretical work, we consider a network of N numbers of variables. An N-dimensional statistical thermodynamics analysis allows fundamental questions about the thermodynamic controls in a coordinated network system to be addressed, including: how to identify the most important components of a regulated network, what is the relative thermodynamic contribution of different network components, and what is the thermodynamic price of homeostasis of a regulated network?

We present an experimentally accessible theoretical framework, which we refer to as Thermo-FP, that demonstrates that the rate of free energy dissipation associated with maintaining a nonequilibrium steady state network of intracellular reactions (i.e. the housekeeping heat) can be determined by the covariances in the temporal fluctuations in the components of the network together with the gradients of the potential landscape. This analysis is enabled by a recent application of FP dynamics [35] to demonstrate that the Boltzmann H-function explicitly connects the relative entropy of a relaxing population with the rate of dissipation of free energy involved in maintaining the network at steady state.

Results

The theoretical framework

In this work, we develop a thermodynamic description of biological networks at the level of microscopic interactions between network variables. Experimentally accessible measurements of network variables at the level of single cells can provide data about the fluctuations in, and dynamic interactions between, those network variables. We build on the rigorous Kullback-Leibler based definition of relative free energy presented by Rao and Esposito [36]. Instead of using a master equation, we apply FP dynamics and the Boltzmann H-function to describe the rate of approach to steady state in terms of the dynamical fluctuating behavior of network components, and the thermodynamic quantities that can be derived from the correlations in fluctuations of the network variables. We call this approach for interpreting the kinetics of phenotype expression Thermo-FP. This coarse-grained approach diverges from a master equation approach, bypassing the need to have detailed knowledge of explicit reaction steps which are often difficult to know with confidence; it allows evaluation of the thermodynamics of complex networks for which there is insufficient knowledge to write chemical rate equations. This approach emphasizes the role of fluctuations or “noise” in controlling a biological system and provides a direct link between the dynamic correlations between network variables and a thermodynamic understanding of network size, composition, and relative importance of network variables as will be shown. This is possible through experimentally measurable quantities.

We present Thermo-FP to describe the evolution to, and maintenance of, a steady state of an N-dimensional network probability distribution with N x N diffusion matrices of variances and covariances in and between network variables. The magnitude of the coupled fluctuations between the different network variables, i.e., the covariances that comprise the diffusion matrix, is a measure of the strength of the physical and functional interaction between those variables [21]. The coordinated relationships between network components create an organized structure, thus reducing system entropy and resulting in the production of dissipative heat [37, 38].

Here we treat multidimensional landscapes of biochemical networks by applying the mathematical properties of positive definite quadratic forms and normal mode analysis common to mathematical physics and statistical mechanics [39]. A diffusion matrix of the dynamic covariances between the multiple variables of the network can be rotated through normal mode analysis to identify complex collective modes of network variables. We show that the product of each eigenvalue of the rotated diffusion matrix with the square of the gradient of the rotated multidimensional landscape allows determination of the contribution each degree of freedom makes to the rate of dissipation of heat that maintains the regulated network of reactants.

We begin by showing that the Boltzmann H function describes the free energy of a population relaxing by Thermo-FP dynamics to a steady state distribution.

The steady state landscape and the Boltzmann H function

The free energy of a system approaching a nonequilibrium steady state can be identified with the Boltzmann H-function, H(t), and is the rate of dissipation of free energy as a system that is transiently perturbed away from its steady state relaxes to a steady state distribution, W.

H(t) is a relative free energy defined by in terms of probability densities, where W1({x},t) is the probability density of the microstates at some time, t during relaxation, W2({x},t+τ) is the probability density at time t+τ (depicted schematically in Fig 2), and T is the thermodynamic temperature of the network system in contact with an isothermal heat reservoir [40]. We assume that no entropy production is associated with temperature gradients. Thus, the relative free energy of the system is defined by the average of the logarithm of the ratio of the occupation probabilities of the microstates, {x}, of the distribution. Keeping in mind that time-dependent relaxation is determined by stochastic fluctuations of characteristic rates, a series of probability density functions is observed over time (Fig 1), each of which can be described by FP dynamics in terms of a potential force, or drift, term and a diffusive term at every microstate; the former corresponds to the gradient of the landscape on which the microstate resides, and the latter is described in an N by N diffusion matrix of variances and covariances of the fluctuations of the N network variables (see S1 Text). As the perturbed system relaxes, H(t) decreases with time () and is a minimum at steady state, SS.

Fig 2

A steady state distribution, if transiently perturbed, will relax back to the steady state distribution.

(Left) The rate of dissipation of free energy during relaxation, , approaches 0 as the steady state (Wss) is reached. A thermodynamically open, nonequilibrium system is less entropic, and of higher free energy, than an equilibrium system. (Right) H(t) corresponds to population distributions (shown here for 1-dimensional distributions). A subpopulation resulting from cell sorting (Day 0) is allowed to relax, and over time, the population achieves the initial steady state distribution from which the subpopulation was taken (3).

The time derivative of the Boltzmann H-function, can be written as a quadratic form [35, 41], consistent with FP dynamics: where is the microscopic relative free energy associated with the microstate {x} of the landscape. The quadratic form in the integrand of Eq 2 is the inner product of the probability current and the thermodynamic force, which explicitly defines the contribution to the rate of free energy dissipation due to fluctuations in network variables. The gradients of the landscape, and are analogous to the gradients of a chemical potential; and D is a diffusion matrix of experimentally determined dynamic covariances between the N variables [41-44]. The N-dimensional diffusion matrix, D, reflects the fluctuations in each of N network components, i.e. the diagonal elements (variances), and the dynamic correlations (covariances) between the N variables, i.e. the off-diagonal elements of the diffusion matrix.

The energetic components of the non-equilibrium steady state

By setting in Eq 2, where W({x}) refers to a time-invariant steady state, we see that Thermo-FP dynamics ensures that vanishes as the stationary state is approached and R approaches unity and ln R ({x},t) approaches 0.

At equilibrium, the rate of increase in entropy and the rate of dissipation of free energy will approach zero. Living biological systems do not exist at equilibrium, and instead reach a nonequilibrium steady state with reduced entropy and higher free energy compared to an equilibrium state. The equilibrium state is defined by detailed balance, in that the probability flux between any pair of points on the landscape is equal to the probability flux in the reverse direction; in the nonequilibrium state there is a lack of reversibility. Regulated interactions between network components create a coordinated set of reactions; this allows the system to respond dynamically to arbitrary perturbations in order to recover homeostasis. This molecular organization reduces entropy and keeps the network system in a nonequilibrium steady state that is maintained by the injection of free energy from outside the network system; that energy is dissipated by the system as heat [37, 38]. During relaxation to an equilibrium steady state, is the rate of free energy dissipated from the production of entropy during approach to the steady state, i.e., where is the entropy production rate. During relaxation to a nonequilibrium steady state distribution, is the rate of free energy dissipated from two sources, namely , and the rate of dissipation of heat from the free energy that is required to maintain the nonequilibrium steady state, , the housekeeping heat. (The dot indicates the time rate of change of the quantity). Thus, for a nonequilibrium system: [38], [38], and dH/dt≤0 [44], and

Before reaching steady state: where is the instantaneous housekeeping heat and is a functional that varies with the probability distribution, W({x},t); W is the nonequilibrium steady state; and W is the equilibrium steady state (where detailed balance holds). The first gradient term in the bracket before the diffusion tensor is the thermodynamic force at steady state and the remainder of the term is the probability flux. The diffusion term, D({x}), explicitly connects network component fluctuations to housekeeping heat.

The instantaneous housekeeping heat as shown in Eq 4 is a bilinear form, which at steady state reaches a constant value:

The steady state expression for is a quadratic form for any trial steady state distribution. The steady state distribution that is a solution of the multidimensional FP equation for the cellular network minimizes the quadratic form in Eq 5.

The Boltzmann H function thus provides us with a measure of the housekeeping heat, introduced by Oono and Paniconi [45] and further developed by others [38, 46]. The positive, semi-definite quadratic form of at steady state is immediately apparent from Eqs 4 and 5.

Normal mode representation of the diffusion matrix

The nonequilibrium steady state distribution of a cell population can be thought of as composed of microstates, {x}, i.e. the different ways the population can achieve its continuum of phenotypic expressions. The phenotypic expressions of the population are determined by the dynamic interactions of the components of the network; these interactions are represented in the N x N diffusion matrix, D({x}), as the autocorrelations and the dynamic covariances of the fluctuations of the activities of the N network variables. The dynamic covariances are the mean square fluctuations about the average relative expression levels divided by time. The H function provides us with the time-derivative of the relative free energy of the system, from which we can identify , the rate at which heat is dissipated during the maintenance of the nonequilibrium steady state of the network. We will use these relationships to determine the relative rate of heat dissipation by the various network components to support the homeostatic steady state distribution.

The normal modes of the diffusion matrix result from rotation of the matrix D({x}) to a diagonal form and provide us with a matrix D({x*}) in which the eigenvalues of all diagonal elements are greater than 0 and the off-diagonal elements are equal to 0, i.e., the components of the matrix are independent of one another. The purpose of this transformation is to identify combinations of network components as composite variables that are the major contributors to the rate of heat dissipation of the system; these are effectively the degrees of freedom. The original diffusion matrix consists of self- and cross-correlations between network variables, where some of the cross-correlations can be positive and some negative. Rotation of the matrix allows us to define composite variables (degrees of freedom) as clusters of N-wise interactions between the N components of the network. The rotation operation guarantees eigenvalues that are positive as well as eigenvectors that are mutually orthogonal, and yields an expression for the energetics of the system at steady state as

In Eq 6, λ({x*}) are the eigenvalues of the rotated diffusion matrix and represent interactions between network components and the magnitude of their coordinated fluctuations. Eigenvectors represents ways, i.e., modes, in which network components are organized with respect to their interactions with one another. The eigenvalues are the diffusion coefficients associated with those independent collective modes of the network components. The term V({x*}) is equal to the gradient of the potential defined by the rotated landscape and expressed as . It should be noted that the eigenvalues as well as the gradients of the potentials are functions of the entire coordinate set of x*. Eq 6 shows that the rate of free energy dissipation to maintain the nonequilibrium steady state will be largest when eigenvalues are large because of strong dynamic interactions between components of the network and when those interactions are occurring in a part of the landscape that is characterized by steep gradients. The rotated diffusion matrix produces positive eigenvalues and corresponding orthogonal eigenvectors, which simplifies the next step which is to sum the most important contributors to the rate of free energy utilization and heat dissipation.

Network contributions to

in Eq 6 is the rate at which the N-dimensional network is dissipating heat associated with the maintenance of the nonequilibrium steady state. It can also be regarded as the rate at which external free energy is injected into the network. The summation term on the right side of Eq 6 is a positive definite quadratic form; in keeping with statistical thermodynamics we use this term to determine the energy that each of the degrees of freedom contribute to the system. The magnitude of λ({x*}) and [V({x*})]2 in integrated over the probability distribution defines the rate of heat dissipation by the various components of the homeostatic network. This treatment implies that some interactions between network components or variables are more dissipative than others and therefore are more important thermodynamic contributors to the steady state. Below we will address the significance of this.

If we consider each composite network variable that is identified by rotation of the matrix as a contributor of a degree of freedom, their sum, C, represents the total rate of heat dissipation associated with the system as shown in Eq 7.

Each term is an implicit function of x* and reflects a heat that is dependent on the magnitudes of the eigenvalues of the matrix, λ, and the square of the gradient, , of the landscape that corresponds to the rotated coordinate system. Together these terms constitute a quadratic representation of the microscopic dissipation of the regulated circuit, a sequence of partial sums that increases monotonically as the number of terms increases. We are proposing that reactions that involve high rates of free energy dissipation are more important to the function of the cell vis a vis stability and adaptability.

We now make an assumption that for a biological system, there will likely be an upper bound, C, on the local rate of dissipative heat produced. It is well established that living systems are sensitive to nonoptimal temperatures [47], so the rate of heat production cannot exceed the rate at which heat can be dissipated without a failure of the system. This upper limit may arise from the temperature associated with heat generation or the limited rate of transport of energy or matter from the environment into the system. The value of the sum of the energetics of the network components thus cannot exceed the upper bound denoted by:

Each of the elements in Eq 8 is a contributor to the summation of network dissipative heat. This treatment thus provides us, in principle, with an experimentally tractable way of assessing the relative contribution that each multivariable component makes to the rate of heat dissipation in the maintenance of the network, and the free energy cost associated with keeping the network in homeostasis. When the number of network dimensions, N, is sufficiently large, C reaches a limit, C and after integrating over the multidimensional probability distribution this constant rate is identical to . A geometric approach for estimating the size of N required for convergence of C to C is described in S1 Text.

The upper limit to the dissipative heat C can constrain the maximum dimension of components in a cellular network. However, at each microstate {x} of the network all the N eigenvalues would not be strictly greater than zero. In fact, if a cellular network is modular then only a small subset of the N network components will have non-zero eigenvalues at each microstate. There is strong evidence that biological networks are hierarchical and modular in their topology [48], which allows for a large number of network components without crossing the dissipation rate threshold.

In addition to the assumption of an upper bound on dissipative heat, we utilize a generalization of the maximum entropy principle to show that the homeostatic heat generation rate, , also has a lower bound [49], i.e., there exists a finite dissipative gap between the non-equilibrium stationary state and the equilibrium (detailed balance) distribution (see S1 Text).

Discussion

It has been frequently noted that stochastic fluctuations in molecular components in individual cells are important to regulatory mechanisms in one-dimensional systems [14, 26, 50, 51] and in multidimensional networks [23, 27, 29–31, 52, 53]. For example, Mojtahedi et. al. [31] analyzed transitions in lineage progression in a population that was attributed to a fluctuation-driven disappearance of an attractor basin. But to our knowledge, the current work is the first to show the direct relationship between fluctuations in and dynamic covariances between network variables and the thermodynamic quantities that contribute to housekeeping heat.

We have limited our theoretical Thermo-FP approach to one that is experimentally tractable. We have shown here how the use of a potential landscape and a diffusion matrix provides a framework for determining the relative energetic contributions of the components of a regulated network. We use FP dynamics because it allows us to apply fluorescence microscopy data from living individual cells and cell populations to directly determine distributions and fluctuations. This approach precludes the need to infer transition probability rates that are needed for a master equation.

With this approach, we derive an experimentally accessible value for , the rate of heat dissipation associated with maintaining the nonequilibrium steady state, directly from the Boltzmann H function. While we have presented this method as applied to analyzing steady state distributions, this approach is also applicable to dynamic population state transitions as discussed in S1 Text.

Large magnitudes of correlated fluctuations of network components are associated with large rates of heat dissipation [54], and our theoretical treatment shows that this is especially true when these fluctuations occur in areas of steep landscape gradients. Noise and complexity are defining features of biological systems, and the Thermo-FP analysis suggests a thermodynamic basis for the relationship between noise and complexity and the stability of a regulated network. Fluctuations are required to ensure ergodicity by allowing cells to escape from deep attractor basins and maintain the stability of the entire landscape structure. We may consider that a biological system requires both stability and adaptability even though these may seem to be opposing characteristics. The Thermo-FP treatment presented here provides a thermodynamic basis for understanding why both the diffusion matrix and barrier gradients are important for maintaining the distribution of phenotypes. Deeper attractor basins associated with non-negligible diffusion coefficients enable stability of the network landscape during nominally constant environmental conditions, and shallower attractors allow adaptation to changing conditions through transition to a new steady state. Very small diffusion coefficients could result in long-lived metastable states, analogous to glassy conditions.

The reasonable assumption that there will be an upper limit to the rate of heat dissipation in a biological system suggests that characteristics of landscapes and diffusion coefficients provide insight into numbers of network variables and the stability and composition of the network. There exists a geometric estimation of the number, N, of variables required for convergence in terms of the volumes of N-balls and N-ellipsoids [55]. Depending on characteristics of the system, such as roughness of the landscape, sufficient numbers of network terms can be predicted to be as small as a few, or many times larger. A reasonable estimate for the number of variables is between 8 and 10. How a sufficient number of terms for convergence to N can be determined is discussed in S1 Text. For example, very large landscape gradients that correspond to large eigenvalues would contribute strongly to the overall dissipative heat of the network, and a small number of such contributors may be sufficient to reach a thermal limit. This condition would be indicated by a landscape containing one or a small number of very deep attractor basins.

The measurement of dynamic correlations of many cellular variables over time in individual cells is in principle achievable with time-resolved fluorescence microscopy of live cells (see S1 Text for details), especially when enabled by automation and advances in handling of large image datasets [56, 57]. Although transcriptomics analysis can probe a larger number of genes compared to live cell imaging, the appropriate interpretation of the relative significance of these changes to network function and their relationships to one another can be ambiguous [21, 33]. Methods like transcriptomics analysis that rely on “snapshots” of populations at single points in time can infer temporal and treatment-dependent relationships between variables, but real time trajectories of changes in gene expression in individual cells, such as is accessible by live cell imaging, can provide unambiguous determination of correlations in stochastic fluctuations between network variables.

In Thermo-FP, stochasticity is captured in the dynamic fluctuations of the network variables. The magnitudes of the correlations in dynamic fluctuations provide a direct measure of the thermodynamics of network component interactions (since the magnitudes of the fluctuations are proportional to kT, Boltzmann’s constant times temperature). Thus the magnitudes of the correlations are proportional to the relative energetic significance of their contribution to the network. The correlations in fluctuations between each pair of N variables of the network are assigned to an N x N diffusion matrix for each microstate. The rate of heat dissipation associated with each composite variable is determined by the magnitude of the fluctuations and the gradient of the landscape. A variable that is a negligible contributor to dissipative heat would be predicted to play a small role in maintaining the network, and could indicate that the variable is coincidental, but not causative, to network function. To achieve unambiguous interpretation, we show here how these putative relationships can be assessed experimentally by directly measuring trajectories in variable space and their covariance.

This analysis has potential practical implications. is a quantitative measure of the rate of heat produced to maintain a regulatory network. The magnitude of can be used to compare the relative thermodynamic cost of different steady-state phenotypic distributions. For example, as a metric it could provide insight into the thermodynamics of different regulatory networks, or the same network functioning in cells from different individuals. It will be a useful metric to guide cell therapy manufacturing conditions, and to guide the engineering of regulatory pathways in synthetic biology applications.

Conclusions

Thermo-FP analysis, through rigorous connection to the Boltzmann H function and the rate of dissipation of free energy of the nonequilibrium steady state, provides a direct relationship between composite network variables and their contribution to the heat of maintaining homeostasis of the network.

The application of the Thermo-FP approach allows dynamical analysis of network interactions and cell states on a continuous landscape. A multidimensional (or multi-variable) landscape that considers the dynamics of network components can provide unique understanding of the correct interpretation of cellular phenotypic indicators in the context of other network components, stable attractor states and rates at which neighboring phenotypic states can be accessed. Furthermore, we have shown how , the dissipative heat required to maintain a multi-variable network can, in principle, be determined from experimentally tractable data consisting of a steady state distribution and a diffusion matrix of dynamic covariances in network variables. Each eigenvalue/eigenvector of the N-dimensional rotated matrix represents a unique and independent cluster of cooperative interactions of network components, and each constitutes a degree of freedom of the network.

Experimental observations of short-time changes in multiple network variables allows determination of the extent to which the time-dependent expression of variables is correlated. The observation of these correlations over time in individual living cells provides confirmation of causative relationships between variables.

This approach to analysis of the multivariable landscape of microstates provides unique insight into the components, paths, and the thermodynamic price associated with maintaining a nonequilibrium regulated network. While we have focused this analysis on the steady state distribution, this approach will also be useful for tracking how cellular populations transition from one steady state to another in response to environmental changes, by helping to identify and quantify interactions between network components during transitions.

a. Fokker Planck dynamics and Thermo-FP framework. b. Experimental considerations. c. Convergence of the homeostatic heat to an upper bound: A geometric interpretation. d. The dissipative heat sustaining homeostasis of a network has a lower bound.

(DOCX)

Click here for additional data file.

29 Oct 2019

PONE-D-19-26315

Properties of a Multidimensional Landscape Model for Determining Cellular Network Thermodynamics

PLOS ONE

Dear Dr. Plant,

Thank you for submitting your manuscript to PLOS ONE. The paper was sent to two reviewers, who have now commented upon its suitability for publication. Both reviewers have raised major concerns that should be addressed before the paper can be considered further. These concerns include, but are not limited to, incomplete reference to previous work, novelty in the light of that previous work, and relevance to biology. Please address all the comments (see below, notice also commented pdf from reviewer 1) in your resubmission.

We would appreciate receiving your revised manuscript by Dec 13 2019 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Jordi Garcia-Ojalvo

Academic Editor

PLOS ONE

Journal Requirements:

1. When submitting your revision, we need you to address these additional requirements. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. We note that your manuscript does not contain a labelled Methods section. Please consider if having a designated Methods section would be beneficial to the reader.

3. We note that Figures 1 and 2 in your submission contain copyrighted images. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.We require you to either (1) present written permission from the copyright holder to publish these figures specifically under the CC BY 4.0 license, or (2) remove the figures from your submission:

1.    You may seek permission from the original copyright holder of Figures 1 and 2 to publish the content specifically under the CC BY 4.0 license. We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:

“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”

Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission.

In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”

2.    If you are unable to obtain permission from the original copyright holder to publish these figures under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: No

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: N/A

Reviewer #2: N/A

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: In their ms, Hubbard and co-workers use an FPE approach and thermodynamic arguments to investigate cell population dynamics and stability. Overall, the ms is written well but some parts would need more clear explanations as stated below and commented in the attached pdf.

While the ms is overall interesting, I have some general concerns with the content and presentation. In particular, their main result for the KL based definition of a relative entropy was derived before and a thermodynamic interpretation of high dimensional readouts was done by others, too. Therefore, the innovation of the ms is a bit limited.

Major points:

- The term Thermo-FP is a bit irritating. The FPE is rooted in statistical physics as an approximation of the Master Eq. and hence per se a TD description (including temperature which was not considered here explicitly). What the authors probably mean with the “Thermo” prefix is a thermodynamic interpretation of the fluctuation correlations? This should be clearly stated and explained in the main text.

- While their argumentation and interpretation is meaningful, the presentation could be improved to make it also readable for non-experts.

- More importantly, I missed some key references and discussions with published work.

o The connection to stochastic thermodynamic (sTD) is not done in the main text (besides one ref to Seifert in another context). In this respect, the statement that “For a nonequilibrium steady state system, the occupation of microstates results from irreversible processes ..” (l. 138) is not correct in sTD processes are reversible!

o The KL based definition of the relative free energy is rigorously derived by Rao et al. (https://journals.aps.org/prx/abstract/10.1103/PhysRevX.6.041064) and should be mentioned and compared here.

o The thermodynamic interpretation of fluctuations wrt to cell states and their stability was also studied by e.g. Chen (https://www.nature.com/articles/srep00342) or Huang (https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.2000640).

- Given these published studies, the innovation of the ms is not really clear (and not clearly stated).

- For me the connection from the physical non-equilibrium state to the biological non-equilibrium state is not clearly described. While from the chemo-physical perspective “ENERGY” would correspond to ATP or similar energy substrates, the biological state is here described by the gene regulatory network state. How are these two levels linked? Is the physical heat dissipation comparable with biological ordering? In this context, Eqs. 3-6 should be explained in more detail to enable plausibility check. Why should Eq. 4 be general valid?

- Given these critical comments above, the proposed application to biological systems is interesting but without a proof of concept application to investigate the network contribution (even with public available data) and a concrete comparison with other methods including those mentioned above, the applicability of the suggested framework cannot be judged.

Besides these major points, I highlighted and commented some minor points in the attached pdf.

Based on this evaluation, I think that the ms is interesting but does not reach the level of PLoSONE in its current form.

Reviewer #2: This manuscript aims to develop a framework for the thermodynamics of a cell population at nonequilibrium steady states. The analysis referred to as Thermo-Fokker-Planck gives insight into the relative contributions of various network components to the relaxation process. The original method was developed for nonequilibrium steady states of nonliving systems. Efforts are made to apply this method to a population of living cells.

These are interesting ideas and the direction is important for today’s physical biology. However, it would be very useful to make the text easier to digest for both non-theoretical physicist and experimental biologist readers unfamiliar with the nonequilibrium thermodynamics bases of the manuscript.

(1) The way this method transfers from nonliving to living systems is unclear. What exactly are the sources and sinks of entropy and energy? In what sense and why are cell populations nonequilibrium? What does temperature mean here? Some clarifications are needed to make all this useful for the community.

(2) While the manuscript follows the spirit of theoretical papers such as Ref. 38 by Oono & Paniconi, it should also incorporate the spirit of Shin-ichi Sasa & Hal Tasaki, Journal of Statistical Physics 125(1), 2006, which should be cited. What Sasa & Tasaki exemplifies is how very simple, realistic systems such as sheared flow or thermal flow can be used to demonstrate the applicability of theory. The same should be done here for at least one or two biological systems: what plays the role of a “wall” (as in sheared flow) for cell populations?

(3) There are some statements that often fail in biological, cellular systems. For example, ergodicity (line 47) and detailed balance (lines 193-195) completely fail if protein levels affect the growth rate (and thereby the dilution rate). That is, the steady-state moments of time courses from tracking single-cell lineages over time will differ from steady-state moments over cell populations at any given time. This effect is described in PMID:22511863 and PMID:30341217, which would be worth citing and discussing. The statements about ergodicity and detailed balance should include the limitation that these are valid only if growth rates do not depend on protein levels.

(4) Line 141: “approaching its nonequilibrium steady state, entropy decreases over time” – this statement should be explained and references should be provided as it is unusual for anyone familiar with standard, classical thermodynamics.

(5) Figure 2: only the high-sorted population is shown over time. In addition, the unsorted population and the low-sorted populations should also be shown at the same time points.

(6) Boltzmann’s constant and temperature do not appear in the formulas of the paper. While it is OK to omit them, their meaning should still be clarified. The approach should be developed with k and T present and then they can be dropped once it is clear what happens with their incorporation. In fact, the temperature here is probably related to the fluctuations of molecule concentrations or cell states, meaning that the temperature may not be identical to the typical “absolute temperature” in statistical physics of nonliving systems. This should be clarified.

(7) It would be helpful if the method could be illustrated on a very simple, 1- or 2-dimensional system, such as a constitutively synthesized protein with or without self-regulation or something similar, using actual matrices, probability distributions, etc.

(8) Related to the previous comment, the heat terms may not be the usual heat measured in nonliving systems. This should be discussed and a practical interpretation for the heat terms should be provided.

(9) The reason for assuming the “upper bound” (line 291) should be clarified. “Temperature associated with heat generation” etc. is unclear because heat and temperature are unclear (see above).

(10) There are some typos throughout the text that should be corrected: “there can significant dynamics variability”,

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: PONE-D-19-26315_reviewer_commented.pdf

Click here for additional data file.

9 Jan 2020

RESPONSES TO REVIEWERS

We would like to begin by thanking the Reviewers for their insightful comments, which have helped greatly to improve the manuscript. We have made a number of changes to the manuscript, in addition to the comments and explanations provided below.

In some cases, we have combined related Reviewer comments to provide a complete and concise response. We have also addressed the comments that Reviewer 1 provided in a pdf mark-up of the manuscript; these changes are embedded within the manuscript.

Responses to Reviewer 1

The term Thermo-FP is a bit irritating…What the authors probably mean with the “Thermo” prefix is a thermodynamic interpretation of the fluctuation correlations? This should be clearly stated and explained in the main text.

Yes, the reviewer is correct, that was exactly what we intended. Thank you for the suggestion. We have made this clear by revising the title of the manuscript and adding additional explanation of the term Thermo-FP in the Abstract and Introduction.

The connection to stochastic thermodynamic (sTD) is not done in the main text (besides one ref to Seifert in another context). In this respect, the statement that “For a nonequilibrium steady state system, the occupation of microstates results from irreversible processes ..” (l. 138) is not correct in sTD processes are reversible! This sentence is ambiguous and has been deleted.

o The KL based definition of the relative free energy is rigorously derived by Rao et al. (https://journals.aps.org/prx/abstract/10.1103/PhysRevX.6.041064) and should be mentioned and compared here.

Our decomposition of dH/dt, the rate of free energy degradation, is conceptually identical to the decomposition performed in previous studies including Rao et al. We decompose dH/dt into adiabatic and non-adiabatic components, namely the housekeeping heat and entropy production rate due to transient dynamics. Our work builds on the rigorous KL based definition of the relative free energy derived by Rao. The work of Rao et al is now cited in Results.

o The thermodynamic interpretation of fluctuations wrt to cell states and their stability was also studied by e.g. Chen (https://www.nature.com/articles/srep00342) or Huang (https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.2000640).

- Given these published studies, the innovation of the ms is not really clear (and not clearly stated).

We apparently have not articulated the innovative aspects of our theoretical framework. These innovative aspects are:

1) Our derivation for the housekeeping heat results in a quadratic form representation that provides insight into the driving force and flux that are the origin of the heat dissipation of the nonequilibrium steady state (NESS). In the form presented in Eq. 4 and 5 (previously Eq. 5 and 6), the driving force for the NESS has a clear interpretation, namely, the difference between the potential of the NESS and the corresponding equilibrium distribution (or detailed balance). The quadratic form also makes it intuitive how the fluctuations in network components, as reflected by the diffusion matrix, contribute to the housekeeping heat of the NESS.

2) The quadratic form also lends itself to normal mode analysis, which allows identification of independent coordinated network components, the dominant modes of relaxation (degrees of freedom) of the network, that contribute to the housekeeping heat (see section “Normal Mode representation of the diffusion matrix”). As a result, this analysis uniquely applies the housekeeping heat to the interpretation of the organization of cellular networks. This implies that there is a minimum heat associated with correlated fluctuations because the quadratic form is bounded by eigenvalues of the diffusion tensor.

3) Fokker-Planck dynamics makes it straightforward to use experimental data quantifying the phenotype of cells and cell populations to analyze complex network analysis. Data on dynamic covariances in phenotype from quantitative live cell microscopy are becoming increasingly available and the corresponding landscape gradients can be obtained from steady-state cell population data.

We have modified the manuscript (Introduction, Results, Discussion and Conclusions) to more clearly highlight these innovative aspects.

We have also integrated the references brought to our attention by the reviewer in the Discussion section.

The work of Chen et al. and Huang et al. and others are in the same spirit as our theoretical approach in that they apply dynamical systems models to interpret transient behavior of cell populations. While it is not unusual to derive thermodynamic quantities from these types of models, our formulation of the housekeeping heat through experimental measurements of dynamic covariances is novel. Furthermore our approach is a method to apply the housekeeping heat to interpret the organization of cell networks. With respect to Chen et al. and Huang et al. our work differs in two important ways:

1) The theoretical approaches of Chen et al and Huang et al use the ensemble averages of cell population distributions, while our approach uses the dynamic covariances in phenotype that can be observed in individual cells, as well as the shape of the population distribution. By measuring changes in phenotype (such as expression of a fluorescent protein) in the same living cell over time, one can determine the fluctuation rate for that phenotypic variable and differences in fluctuation rates in different cells. By measuring more than one variable in the same cell, one can measure the correlations in fluctuations for those variables in individual cells, and the differences in correlation in different cells that may be expressing different amounts of one or the other of the variable. Access to dynamic covariance through quantitative live cell microscopy is becoming increasingly available and theoretical frameworks that can be used to interpret cell signaling networks are needed.

2) We show how the shape of the distribution and the dynamic covariances are related to the housekeeping heat required to maintain the NESS (see description of innovative aspects above). Neither Chen et al. nor Huang et al. make connections to this thermodynamic quantity.

- For me the connection from the physical non-equilibrium state to the biological non-equilibrium state is not clearly described. While from the chemo-physical perspective “ENERGY” would correspond to ATP or similar energy substrates, the biological state is here described by the gene regulatory network state. How are these two levels linked? Is the physical heat dissipation comparable with biological ordering? In this context, Eqs. 3-6 should be explained in more detail to enable plausibility check. Why should Eq. 4 be general valid?

Our mathematical framework is applicable to both physical and biological nonequilibrium states. The constraints on the biological system could be considered to be equivalent to boundary conditions.

The dissipative heat is indeed a result of biological ordering and the reduction of entropy. This concept is well developed by Prigogine (citation included in the text). The network is maintained by injection of energy sources or other biochemicals into the system to drive the reactions that maintain the organization of network components, and dissipative heat is released. The advantage of coarse graining with FP is that the molecular species that are involved in maintaining the network variables at their appropriate levels don’t have to be stipulated. In fact, the biological data that are collected (for example by fluorescence microscopy) provide measurable quantification of the response of network variables that could be the result of many upstream (unidentified) factors.

The levels of network variables are measured directly in each individual cell over time; no assumptions of how the microstates are reached are necessary. The coordinated fluctuations between network variables provide the thermodynamic quantities associated with their interdependence as a function of microstate on the landscape. Eq. 4 and 5 (formerly 5 and 6) show how the coordinated fluctuations in network variables (Dij) contribute to the housekeeping heat required for maintaining the nonequilibrium steady state of the network.

In order to connect more seamlessly with other authors, We have removed the term Q.ENERGY and replaced it with Q.HK, the dissipative heat that remains associated with a nonequilibrium steady state after (dH(t))/dt has reached zero. We also modified Eq. 3 so that it is directly comparable to the housekeeping heat described by previous authors (Oono and Paniconi).

- Given these critical comments above, the proposed application to biological systems is interesting but without a proof of concept application to investigate the network contribution (even with public available data) and a concrete comparison with other methods including those mentioned above, the applicability of the suggested framework cannot be judged.

In Sisan et al (PNAS 2012), we showed how experimental data can be used to compute gradients of the landscape and diffusion coefficients, and therefore, in the context of this work, how to compute housekeeping heat. This work follows on that study (which was on a 1-dimensional system) and describes how experimental data from multiple network variables can be collected and interpreted with respect to thermodynamic quantities. In order to make this connection between this work and the previous work, we have elaborated on the methods and results of Sisan et al in the Introduction. In that work, we showed that subpopulations that were each expressing different levels of the reporter relaxed to the steady state landscape with very different kinetics. Langevin/FP dynamics and the diffusion coefficient calculated from experimentally measured fluctuations in fluorescent protein expression allowed nearly perfect prediction of the relaxation kinetics in the absence of any adjustable parameters. This fairly remarkable finding encouraged us to develop the Thermo-FP framework for a multi-variable network system.

Responses to Reviewer 2

This manuscript aims to develop a framework for the thermodynamics of a cell population at nonequilibrium steady states. The analysis referred to as Thermo-Fokker-Planck gives insight into the relative contributions of various network components to the relaxation process. The original method was developed for nonequilibrium steady states of nonliving systems. Efforts are made to apply this method to a population of living cells.

The current theoretical work is based on previous experimental work (Sisan et al, 2012 PNAS) of the cells in a living population. Thus, all our efforts are applicable to living systems. In order to make the connection between this work and previous work, we have elaborated on the methods and results of Sisan et al in the Introduction.

(1) The way this method transfers from nonliving to living systems is unclear. What exactly are the sources and sinks of entropy and energy? In what sense and why are cell populations nonequilibrium? What does temperature mean here? Some clarifications are needed to make all this useful for the community.

The Results section on The steady state landscape and the Boltzmann H function has been substantially edited to focus better on the thermodynamic aspects of the nonequilibrium steady state that is of primary importance here.

In this work, sources and sinks are the chemical and biochemical conditions that break detailed balance and are the origin of housekeeping heat. In order for cellular networks to maintaining homeostasis, sources of energy or other biochemicals are injected into the system and dissipative heat is released to the exterior of the system as a result. We can speculate about the identity of sources and sinks, but FP coarse graining allows analysis of the thermodynamic contributions of the network without knowing those details. The dynamic data that are collected from living cells over time (for example by fluorescence microscopy) provide a quantification of the response of multiple network variables simultaneously; these responses could be the result of many upstream (unidentified) factors. The advantage of coarse graining with FP is that the molecular species that are involved in maintaining the network variables at their appropriate levels (which are largely unknown) don’t have to be stipulated to gain predictive insight into the parameters of the network.

We have added a citation in the introduction to support the statement: Living cells are a clear example of a nonequilibrium system (Wang 2015 Adv Phys). This citation provides many examples of nonequilibrium thermodynamics studies of cells and cell populations.

In addition, we now make clear that by temperature we mean a real thermodynamic temperature; T is now explicitly written into the equations throughout the manuscript. We assume our systems are isothermal and that no entropy production is associated with temperature gradients. We have included this in the paper in Results.

(2) While the manuscript follows the spirit of theoretical papers such as Ref. 38 by Oono & Paniconi, it should also incorporate the spirit of Shin-ichi Sasa & Hal Tasaki, Journal of Statistical Physics 125(1), 2006, which should be cited. What Sasa & Tasaki exemplifies is how very simple, realistic systems such as sheared flow or thermal flow can be used to demonstrate the applicability of theory. The same should be done here for at least one or two biological systems: what plays the role of a “wall” (as in sheared flow) for cell populations

The concentration of biomolecules that break detailed balance (i.e. the energy sources and starting materials) are analogous to the “wall” in a sheared flow system. The constrained concentrations of biochemical species constitute the sources and sinks in our system.

Regarding an example for demonstration, we refer the Reviewer to Sisan et al (PNAS 2012). In that work, we showed how experimental data can be used to compute gradients of the landscape and diffusion coefficients, and therefore, in the context of this work, how to compute housekeeping heat. This work follows on that study (which was on a 1-dimensional system) and describes how experimental data from multiple network variables can be collected and interpreted with respect to thermodynamic quantities. In order to make this connection between this work and the previous work, we have elaborated on the methods and results of Sisan et al in the Introduction.

(3) There are some statements that often fail in biological, cellular systems. For example, ergodicity (line 47) and detailed balance (lines 193-195) completely fail if protein levels affect the growth rate (and thereby the dilution rate). That is, the steady-state moments of time courses from tracking single-cell lineages over time will differ from steady-state moments over cell populations at any given time. This effect is described in PMID:22511863 and PMID:30341217, which would be worth citing and discussing. The statements about ergodicity and detailed balance should include the limitation that these are valid only if growth rates do not depend on protein levels.

We agree with the reviewer that this is an important point regarding the application of some statistical physics concepts (i.e. ergodicity or detailed balance) to cell populations. Therefore, we have added the following paragraph to the Supporting Information 1B :

“Individual cells in a steady-state population are simultaneously transitioning between microstates and dividing. When the division rate depends on the microstate of the cell (e.g. certain phenotypes have different division rates than others), the moments of the steady-state distribution, W(x), will differ from the moments derived by following a number of individual cells over time (Ref PMID:22511863 and PMID:30341217). A microstate-dependent division rate also means that the landscape derived from the observed steady-state distribution, ln (W(x)), includes the effect of microstate-dependent division rates. An approach for compensating the landscape shape to account for microstate dependent division rates for modeling purposes was shown by Sisan et al (ref Sisan, PNAS, 2012). The landscape transformation was possible because the specific form for the dependence of the division rate on the microstate was identified by quantitative live cell microscopy.“

(4) Line 141: “approaching its nonequilibrium steady state, entropy decreases over time” – this statement should be explained and references should be provided as it is unusual for anyone familiar with standard, classical thermodynamics.

Thank you for catching this typographical error. This sentence now reads “approaching its nonequilibrium steady state, entropy increases over time”.

(5) Figure 2: only the high-sorted population is shown over time. In addition, the unsorted population and the low-sorted populations should also be shown at the same time points.

(7) It would be helpful if the method could be illustrated on a very simple, 1- or 2-dimensional system, such as a constitutively synthesized protein with or without self-regulation or something similar, using actual matrices, probability distributions, etc.

Our previously published experimental work (Sisan et al PNAS 2012) which inspired this current theoretical work is now discussed in greater detail in the Introduction. The full data set with unsorted and the low-sorted populations is presented in that paper which is referenced. That work provides the best example of this approach. The experimental data was collected with fluorescence microscopy and flow cytometry, and the time-dependent analysis on live cells provided the data from which fluctuation rates were calculated. The paper demonstrates how the data provided both probability distributions and diffusion coefficients. These parameters alone provided a remarkably accurate prediction of the rate at which subpopulations of cells relaxed to the steady state. That study was of a one-variable system, and this theoretical work demonstrates how to build on that with a multidimensional network of associated variables.

(6) Boltzmann’s constant and temperature do not appear in the formulas of the paper. While it is OK to omit them, their meaning should still be clarified. The approach should be developed with k and T present and then they can be dropped once it is clear what happens with their incorporation. In fact, the temperature here is probably related to the fluctuations of molecule concentrations or cell states, meaning that the temperature may not be identical to the typical “absolute temperature” in statistical physics of nonliving systems. This should be clarified.

(8) Related to the previous comment, the heat terms may not be the usual heat measured in nonliving systems. This should be discussed and a practical interpretation for the heat terms should be provided.

Temperature and kB are now explicitly written into the equations throughout the manuscript. We assume that our systems are isothermal with a reservoir and that no entropy production is associated with temperature gradients. We have included this in the paper after Eq. 1. The heat referred to is the usual heat measured in nonliving systems. The expectation is that this will be small and difficult to measure explicitly but is identical to the thermodynamic quantity.

(9) The reason for assuming the “upper bound” (line 291) should be clarified. “Temperature associated with heat generation” etc. is unclear because heat and temperature are unclear (see above).

From a practical point of view it is reasonable to assume that there will be a limit to the rate of heat dissipation. While no physical system can withstand an infinite thermal gradient, it is well established that biological systems can be strongly influenced by temperature (Charlebois et al PNAS 2018). This citation has been added in Results.

Submitted filename: Response2Reviewers.docx

Click here for additional data file.

29 Jan 2020

PONE-D-19-26315R1

The Role of Fluctuations in Determining Cellular Network Thermodynamics

PLOS ONE

Dear Dr. Plant,

Thank you for submitting your manuscript to PLOS ONE. As you'll see from the report below, reviewer 1 still has a minor issue with the terminology that I would ask you to address before final acceptance. Upon resubmission I will send the manuscript again to reviewer 1 for a quick final decision.

We would appreciate receiving your revised manuscript by Mar 14 2020 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Jordi Garcia-Ojalvo

Academic Editor

PLOS ONE

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: (No Response)

Reviewer #2: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: N/A

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: While the authors have addressed most of my comments in a satisfactory manner and the ms has improved, I still do not agree with their general statement to reviewer 2 that "In addition, we now make clear that by temperature we mean a real thermodynamic temperature; T is now explicitly written into the equations throughout the manuscript." In best case, their temperature is a relative temperature that can be compared between the different system states but cannot be an absolute temperature (measured in K) because the underlying molecular mechanisms are not resolved. Hence, investigating 2 different systems or 2 different measurement wrt monitored proteins by microscopy of the same system will/can lead to different "temperatures" that can only barely be compared. This should be clarified in the text to avoid confusion of the readership. Otherwise, the ms is in a solid shape now.

Reviewer #2: I would like to thank the Authors for addressing my questions. I would like to recommend publication.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

30 Jan 2020

We have made an additional edit to the manuscript to address the concern of Reviewer 2. We have made a modification to our definition of T in Eq 1 and cited this use by Jarzynski in his seminal 1997 work.

“…and T is the thermodynamic temperature of the network system in contact with an isothermal heat reservoir (40). We also assume that no entropy production is associated with temperature gradients. “

We trust that making it clear that the network is in contact with a thermal reservoir at equilibrium is the point the reviewer wanted us to make.

Submitted filename: Response to Reviewer.2.docx

Click here for additional data file.

21 Feb 2020

The Role of Fluctuations in Determining Cellular Network Thermodynamics

PONE-D-19-26315R2

Dear Dr. Plant,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Jordi Garcia-Ojalvo

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: (No Response)

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: (No Response)

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: (No Response)

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: (No Response)

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

28 Feb 2020

PONE-D-19-26315R2

The Role of Fluctuations in Determining Cellular Network Thermodynamics

Dear Dr. Plant:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Jordi Garcia-Ojalvo

Academic Editor

PLOS ONE