Cellular Compartmentalization as a Physical Regulatory Mechanism of Signaling Pathways

Abstract

Cells compartmentalize biochemical processes using physical barriers in the form of membranes. Eukaryotes have a wide diversity of membrane-based compartments that can be used in this context, with the main ones being the extracellular membrane, which separates the inside from the outside of the cell, and the nuclear membrane, which separates the nucleus from the cytoplasm. The nuclear membrane not only isolates and protects the DNA and the transcription and replication processes from the other processes that are occurring in the cytoplasm but also has an active role in the regulation of cellular signaling. The TGF-$\beta$ pathway is one of the most important and conserved signaling cascades, and it achieves compartmentalization using a well-tuned balance between the import and export rates of the active and inactive forms of key proteins. Thus, compartmentalization serves as an additional regulatory mechanism, physically isolating transcription factors from their targets, influencing the dynamics and strength of signal transduction. This contribution focuses on this biophysical layer of regulation, using the TGF-$\beta$ pathway to illustrate the molecular mechanisms underlying this process, as well as the biological consequences of this compartmentalization. We also introduce a simplified mathematical formulation for studying the dynamics of this process using a generalized approach.

1. Introduction

The development and organization of multicellular organisms are mainly regulated by intercellular communication in the form of signaling molecules [1]. These molecules trigger intracellular events that are processed by the cell in the form of signaling pathways, which ultimately lead to changes in the expression, regulation, and dynamics of proteins. Most of these regulatory events are biochemical in nature, e.g., the reversible covalent addition of a phosphate group (phosphorylation), sugar molecules (glycosylation), ubiquitin proteins (ubiquitination), acetyl groups (acetylation), and methyl groups (methylation). These post-translational modifications can alter many aspects of the target proteins: phosphorylation can activate or deactivate enzymes and signaling proteins, acetylation of histones can alter DNA–protein interactions, glycosylation can affect protein folding and stability, methylation can affect the dynamics of protein–protein interactions, and ubiquitination can target a protein for degradation [2].

In addition to these biochemical mechanisms, there are other types of regulatory mechanism that involve physical constraints. One of the most studied is protein compartmentalization, i.e., intracellular proteins can be directed to specific organelles or regions within the cell, such as the nucleus or the cell membrane. In addition to providing and maintaining organization at the cellular level, this physical separation of processes is highly dynamic (proteins can change location due to specific signals) and has a strong biological regulatory impact [3]. For instance, constraining a given enzyme to a specific region results in a localized increase in concentration, enhancing its efficiency. One of the most important cases is the compartmentalization of glucokinase, which can shift location in response to substrate availability [4]. In brief, a fraction of the glucokinase can translocate from the nucleus to the cytoplasm after stimulation, resulting in a reduction in its nuclear to cytoplasmic ratio. Another important protein subjected to compartmentalization is Rubisco, a crucial enzyme in photosynthesis and carbon fixation. The compartmentalization of this protein in high concentrations inside the chloroplasts of bundle-sheath cells or in the carboxysomes of cyanobacteria has been shown to affect its catalytic efficiency and, consequently, the overall photosynthetic rate of plants [5]. Interestingly, compartmentalization may restrict the conformational variability in a protein, shaping the evolution of all the living organisms that harvest energy via photosynthesis [6].

Some proteins are known to change function depending on the compartment to which they are confined. For instance, ERp57, a protein that belongs the protein disulfide isomerase (PDI) family [7] regulates protein folding while in the cytoplasm, [8,9]. When the protein is inside the nucleus, its main function is as a transcription factor [10,11,12]. This switch illustrates the importance of localization and its impact in function, and it probably affects many other proteins [13,14,15,16].

In our laboratory, we have focused on this concept of compartmentalization in the context of dimensionality reduction [17,18,19,20]. This concept refers to the transition from moving in a three-dimensional space (the cytoplasm) to a two-dimensional space where molecules are embedded or associated with the cell membrane [20]. This shift from 3D to 2D movement results in a important change in local concentration, increasing the probability of molecular encounters and interactions by orders of magnitude and, consequently, the rate of biochemical reactions [19].

This way, the association of a protein or ligand to a receptor in the membrane catalyzes subsequent reactions by restricting the movement of the ligand to the membrane surface, resulting in an effective increase in its local concentration because the volume in which the molecules can move is significantly reduced; this mechanism of action drives the dynamics of activation of erythropoietin (EPO) [20]. It is the basis of our design of chimeric drugs that can selectively affect tumor cells versus normal cells [17,18,19].

The bidirectional movement of transcription factors between the nucleus and cytoplasm in eukaryotic cells is a common example of intracellular compartmentalization. This exchange of proteins that are dynamically excluded from the cytoplasm or the nucleus is essential for various cellular functions, including gene expression, signal transduction, and cell cycle regulation. The dysregulation of nuclear translocation has been shown to be directly involved in a variety of pathological conditions, with a critical role in maintaining cellular homeostasis.

One of the most important pathways that take advantage of this nucleocytoplasmic shuttling to modulate signal transduction is transforming growth factor $\beta$ (TGF-$\beta$), involved in the regulation of a wide range of signaling functions, such as tissue-specific control of development, morphogenesis, cell proliferation and differentiation, tissue homeostasis, regeneration, and even cell-specific or tissue-specific motility [21]. It has also been shown to have a growth-inhibitory effect on epithelial cells, suggesting a tumor suppressor role in carcinomas and cell fate determination [22].

In this review, we focus on the effect of transcription factor compartmentalization in the context of the TGF-$\beta$ pathway and how this physical regulation is utilized to modulate signal transduction. We start with a review of the main aspects of the pathway, followed by several key studies that investigated how this nucleocytoplasmic translocation shapes the dynamics of the signal. We finalize with an illustration of how mathematical modeling can be applied in the context of protein compartmentalization to study these types of biophysical aspects of biological regulation.

2. The TGF-$\beta$ Pathway Signaling Dynamics and Regulation

The TGF-$\beta$ signaling pathway takes its name from the main extracellular signal that activates the cascade of downstream events: the human transforming growth factor $\beta$ (TGF-$\beta$) ligand superfamily. These dimeric polypeptide growth factors involve thirty-three secreted cytokines [23,24], and the different family members have been given various names according to their molecular identification history: TGF-$\beta$s (three isoforms of TGF-$\beta$s: TGF-$\beta$1, TGF-$\beta$2 and TGF-$\beta$3), activins, bone morphogenetic proteins (BMPs), growth differentiation factors (GDFs), nodal, and Müllerian inhibiting substance (MIS) [25]. Interestingly, this superfamily of signaling molecules is one of the most conserved throughout the animal kingdom, and it is believed to have originated in the early stages of multicellular evolution [26].

Canonical TGF-$\beta$ signaling starts with the reversible binding of one of these dimeric ligands to two transmembrane receptors (type-I receptor or TGF$\beta$RI and type-II receptor or TGF$\beta$RII). After binding, these receptor molecules undergo a conformational change that facilitates their interaction with two other transmembrane receptors (TGF$\beta$RI) [1]. The result is the formation of a heterotetrameric complex composed of two TGF$\beta$RII and two TGF$\beta$RI molecules, represented as (TGF$\beta$RII)–(TGF$\beta$RI) [27]. Once the receptor complex is formed [28], the type-II receptors phosphorylate and activate the type-I kinase receptors, following endocytosis inside vesicles characterized by the presence of the Smad anchor for receptor activation protein (SARA) on their membrane. Once inside the endosomes, the type-I receptors are able to phosphorylate and activate downstream proteins (Figure 1).

The Smad proteins that give name to the SARA protein are themselves the main effectors of the signaling events mediated by TGF-$\beta$ ligand stimulation. The Smads are a family of transcription factors (i.e., they bind directly to DNA) highly conserved throughout evolution, and they come in three types. The first type includes R-Smad (receptor-regulated Smads), which includes Smad1, 2, 3, 5, and 8 [26]. The second type is Co-Smad (common mediator Smads) or Smad4, which serves as a cofactor of the R-Smads [30,31]. The third type is I-Smads (inhibitory Smads), including Smad6 and 7, which inhibit R-Smads’ activation by interfering with their phosphorylation by the receptor complexes [30,32].

3. R-Smad Nucleocytoplasmatic Shuttling as a Modulation of Signal Transduction

As mentioned before, the R-Smads constitute the direct substrates of TGF-$\beta$ family receptors and are the key mediators of the TGF-$\beta$ pathway, translating receptor activation into changes in gene expression levels. All R-Smads share a common molecular structure, consisting of two conserved globular domains: N-terminal MH1 (Mad Homology 1), mainly responsible for direct DNA binding; and C-terminal MH2 (Mad Homology 2) domain, which is mainly responsible for receptor interaction. Both domains are connected by a linker region that is not conserved but contains multiple phosphorylation sites for various kinases that regulate the function and activation of the different R-Smads. Two main sets of R-Smads can be identified depending on their specific receptor: Smad2 and Smad3 are activated via TGF-$\beta$ and activin/nodal [33], while the rest are mainly activated by BMP receptors [26,29,34,35,36].

One of the main characteristics of R-Smads is their ability to translocate between the cytoplasm to the nucleus, depending on their phosphorylation state and organization as monomers or complexes. The sequence of events in this well-characterized process, commonly referred to as Smad nucleocytoplasmic shuttling, occurs as follows: (a) unphosphorylated R-Smads are mainly located in the cytoplasm as monomers; (b) after activation by the type-I receptor, two R-Smads typically assemble with one Co-Smad in the form of heterotrimers, also known as the activated Smad complex (the complex is stabilized by the phosphorylated serines at the C-terminal); (c) these complexes then translocate into the nucleus where they act as transcription factors [37,38]; (d) these complexes are then dephosphorylated via interaction with PPM1A and other phosphatases [39], and this dephosphorylation has been linked to the destabilization of the trimer complex complex [40]; (e) once dephosphorylated, the complexes disassemble; (f) finally, the monomeric R-Smads and Smad4 are exported separately back to the cytoplasm by distinct mechanisms [40].

The two main processes involved in this shuttling are nuclear import and nuclear export through the nuclear pore complex (NPC), a large protein structure embedded in the nuclear envelope that serves as a gateway for the selective transport of molecules between the nucleus and cytoplasm. In brief, the NPCs act as hydrophobic channels [41,42] that allow proteins smaller than 20–30 kDa to pass via diffusion, while larger proteins require active transport mediated by importin and exportin receptors [43,44,45].

4. Dynamics and Regulation of Export and Import of the R-Smads

In recent years, several studies have used different levels of R-Smad protein in the nucleus or cytoplasm to study the mechanisms that regulate this process [39]. These studies have shown that the import activity of the R-Smad complexes relies on the interaction of importins with the MH1 domain (specifically importin-$\beta$1 [46] but also importin 7 and 8 [47]) in a process that depends on the RanGTP/GDP gradient across the nuclear membrane. On the other hand, it has been proposed that Smad2 can also enter the nuclei via the direct interaction of proteins of the NPCs with its MH2 domain [40]. Additionally, several experimental studies have propose that not only multimeric Smads (i.e., as part of complexes) but also monomeric Smads can enter the nuclei t through a mechanism different from that of the R-Smad complexes [40].

Regarding nuclear export, it has been proposed that R-Smad monomers exit the nucleus mediated by exportin 4 [26,48,49] but not exporting 1 (the target of leptomycin B [50]). The export mechanism of the R-Smads has been shown to depend on a Ran-GTPase and Ran-BP3, which recognize unphosphorylated R-Smad protein, resulting from the activity of nuclear R-Smad phosphatases, such as PPM1A [51]. Similar to the import process, Smad2 has been proposed to exit the nucleus via the direct interaction with the proteins of the NPCs [52,53]. This is consistent with the lack of the Lys-rich KKLKK sequence in its MH1 domain, which is critical for the nuclear export of Smad3 [54].

As a result of this complex regulation, the R-Smads in the cytoplasm are mainly in the form of unphosphorylated monomers, while, after activation, they primarily locate inside the nucleus, mostly in the form of complexes. The concentration of R-Smas monomers or complexes in each compartment is determined by the differences in their import and export rates. This way, a higher concentration in the cytoplasm of unphosphorylated R-Smad monomers results from a slower import rate and a faster export rate, compared to the complexes [55]. Each of these rates depends on the ability of R-Smads to interact with NPCs, importins, and exportins in monomer or complex configuration [55,56].

On the other hand, the Smad4 monomer can be found in both the nucleus and the cytoplasm [57] due to its similar rates of movement into and out of the nucleus, while, when part of a complex, it is mainly located inside the nucleus, following the dynamics dictated by the other components of the complex (R-Smads) [58]. Interestingly, other studies showed that the inhibition of CRM1 or exportin 1 leads to the rapid accumulation of Smad4 in the nucleus, suggesting that Smad4 is also actively translocating [50]; therefore, the mechanisms of Smad4 translocation and its impact in the translocation of the active complex needs further study.

5. Taking Advantage of Nucleoplasmatic Shuttling to Study the TGF-$\beta$ Pathway

An advantage of these differential R-Smad levels in cellular compartments is that they can be used as simplified readouts of changes in TGF-$\beta$ pathway activation, measurable by microscopy. These type of studies rely on the direct visualization of the R-Smads or Co-Smads with fluorescence or photoactivatable fusion proteins [59] or the use of highly reliable and specific antibodies that can be used to measure changes in relative levels and compare the experimental conditions in fixed samples. For instance, the seminal contribution by the laboratory of Caroline Hill [56] used FLAG-tagged Smad4 and antibodies against Smad2 and Smad3 to show that the inhibition of CRM1-mediated nuclear export results in very fast Smad4 nuclear accumulation, while endogenous Smad2 and Smad3 remain completely unaffected. The same laboratory [39,60] developed a fusion of Smad2 with GFP to follow the difference in translocation in unstimulated versus TGF-$\beta$-stimulated cells in vivo. This study showed that the nuclear accumulation of the active form of Smad2 is due to a decrease in the rate of the nuclear export of the protein rather than an increase in import [39]. Moreover, it suggested that the import rates of phosphorylated R-Smads (in complex configuration) and unphosphorylated R-Smads (mainly in monomer configuration) are similar, indicating that import proteins act similarly in the activated and deactivated states. If this is correct, the import rate of R-Smads as complexes should be twice as fast as that of the R-Smads as monomers, since each complex is composed of two molecules of R-Smad. On the other hand, the fact that the rate of export of the phoshorylated form of R-Smads is much smaller than the rate of export when unphosphorylated suggests that the configuration as complexes prevents binding with exportin 4.

Another important example that utilizes nucleocytoplasmic shuttling was reported by Liu et al. [61], where they studied how variable dose levels of ligands are translated into intracellular signaling dynamics and how continuous ligand doses can be translated into discontinuous cellular fate decisions. In this study, they found that the pathway processes continuous and pulsating TGF-$\beta$ ligand levels differently, displaying varied sensitivities to short and long pulses, also finding that it displays different sensitivities in response to short and long pulses of ligand. These two types of dynamic responses may be crucial for the role of the TGF-$\beta$ pathway in mediating stem cell differentiation.

The differential concentration of R-Smads as a readout was also used [62] to study the temporal signal by the TGF-$\beta$ pathway and its implications for embryonic patterning. This time, the authors took advantage of the nuclear accumulation dynamics of a GFP-Smad4 fusion protein to measure TGF-$\beta$ pathway activity in response to different types of ligand stimulation. The same fusion protein was used to show pulsatile behaviors in Smad4 nuclear localization [58]. Similarly, quantitative measurements of Smad2/3 nuclear localization illustrated that the TGF-$\beta$ pathway exhibits fold-change detection properties [63], i.e., cells respond to relative rather than absolute levels of TGF-$\beta$. In another context, the accumulation of R-Smads in the nucleus was the main readout used to study the crosstalk between oncogenic Ras and TGF-$\beta$ signaling [64], since Ras activation inhibits the TGF-$\beta$-induced nuclear accumulation of Smad2 and Smad3.

In our laboratory, we have studied how different R-Smad accumulations influence the role of the TGF-$\beta$ signaling cascade during vertebrate spinal cord formation [65,66]. In these studies, we quantified the changes in the nuclear localization of the two R-Smads to explore the effects of the interplay between Smad2 and Smad3. We showed that the heterotrimers Smad2–Smad3–Smad4 induce cooperation and antagonism in the activation of the targets of the two R-Smads and that this dual role strongly influences the balance between neural stem cell proliferation and differentiation [65]. Using immunofluorescence and confocal imaging, we visualized and quantified R-Smad expression patterns at different developmental stages, showing pathway activation (i.e., R-Smads become more nuclear) in the transition zone (where stem cells cycle and perform interkinetic nuclear migration) and the differentiated zone (where differentiated neurons migrate after exiting the cell cycle). The study was complemented by numerical simulations that allowed us to interpret the dual antagonism and cooperation as a consequence of Smad heterotrimer assembly [65].

6. Mathematical Modeling of Regulatory Compartmentalization

Mathematical modeling in molecular biology is a complementary approach to the more time-consuming and expensive experimental approach, which can help explain results, make predictions, and test hypotheses that complement in vivo and in vitro studies. In the context of TGF-$\beta$ signaling, several key contributions have combined modeling approaches with experimental results to try to understand the different aspects of the regulation of compartmentalization as a physical regulation of gene expression [61,65,67].

One of the earliest approaches to studying the TGF-$\beta$ pathway dynamics using numerical simulations focused on the dynamics in experiments of fluorescent recovery after photobleaching (FRAP) [67]. The results support the hypothesis that there is a dynamic maintenance of Smad nuclear accumulation during signaling [67]. In another study, the authors developed a model that shows how cells can integrate repeated pulses and maintain the phosphorylated Smad2 concentration [61]. The main conclusion derived from the model is that the pathway is able to respond with different sensitivities to different ligand doses and time scales [61].

One main caveat of these early approaches is that the models used are composed of around 20 ordinary differential equations, making them very difficult to interpret. More recent examples have opted for more simplified but less general mathematical formulations, for instance, using the “pseudo-steady-state” approximation or “rapid equilibrium assumption” to reduce the network signaling reactions based on an analysis of the different time scales of the individual reactions [68]. This model predicted a critical role of potential positive and delayed negative feedback loops on the regulation of the TGF-$\beta$ signaling system.

In line with simplified models, composed of the minimal set of equations required for explaining experimental observations, we developed a mathematical model based on differential equations to support the hypothesis that the heterotrimer Smad2–Smad3–Smad4 is formed during vertebrate neurogenesis and that its formation could explain the fact that Smad2 and Smad3 antagonize in some targets while cooperate in activating the transcription of others [65]. The model prediction was then validated using immunoprecipitation [65].

In this section, we propose, using the same philosophy of simple modeling approaches, a very simple yet powerful mathematical approach to study the dynamics of this nucleocytoplasmic shuttling. The model was designed with ordinary differential equations derived using the mass action law and was used to reproduce the dynamics of a R-Smad protein in and out of the nuclei. Using this simplified system, it was easy to study the main mechanisms and parameters that affect this dynamic equilibrium without needing to perform complicated calculations or simulations.

To build the model, let us consider an R-Smad protein inside of a cell that can enter and exit the nucleus. These two potential states of the R-Smad protein are represented in the model as C (for cytoplasmic) and N (for nuclear). Based on the molecular details of the process reviewed above, we know that the processes crossing the NPCs are actively mediated by interactions with other proteins, which we call $Imp$ (for import) and $Exp$ (for export). This way, to enter into the nucleus, a cytoplasmic molecule C must be linked to the import enzyme $Imp$ and forms a complex $C|Imp$. This complex is then allowed to enter the nucleus. Once in the nucleus, the $C|Imp$ complex is disassembled into one molecule N (since the R-Smad molecule is now inside the nucleus) and the $Imp$ enzyme.

Similarly, the export process follows the same rationale: The nuclear molecule N binds to the export enzyme $Exp$, forming the $N|Exp$ complex. This complex is now able to cross the NPCs and translocate to the cytoplasm. Once in the cytoplasm, the complex is disassembled, giving one C (our molecule is now in cytoplasmic configuration again) and the enzyme $Exp$. The rates of the nuclear import and export processes for $C|Imp$ and $N|Exp$ complexes are $k_{imp}$ and $k_{exp}$, respectively. The formation rates of $C|Imp$ and $N|Exp$ are $k_1$ and $k_2$, respectively. We included the fact that complex formation is a reversible process in the model, which means that the complex can spontaneously disassemble before crossing the NPCs. The dissociation rates of $C|Imp$ and $N|Exp$ are labeled as $k_{-1}$ and $k_{-2}$, respectively. Moreover, the dissociation of each complex once it crosses the NPCs is considered to be instantaneous, and $Imp$ and $Exp$ transport enzymes that can translocate freely between the nucleus and cytoplasm (they are not limiting reagents) [69,70]. The whole system of interactions can be illustrated using a chemical reaction formalism as

$\begin{array}{cc}\hfill \mathrm{Import}\text{:}& C+Imp\underset{k_{-1}}{\stackrel{k_1}{\rightleftharpoons }}C|Imp\stackrel{k_{imp}}{\to }N+Imp\\ \hfill \mathrm{Export}\text{:}& N+Exp\underset{k_{-2}}{\stackrel{k_2}{\rightleftharpoons }}N|Exp\stackrel{k_{exp}}{\to }C+Exp\end{array}$

The scheme of the interaction is illustrated in Figure 2.

6.1. Using the Mass Action Law to Model a Biological System

To translate this system of interactions into differential equations, we use the mass action law formalism, a fundamental principle in the mathematical modeling of chemical and biological systems. In brief, the law of mass actions states that the rate of change of a species involved in a chemical interaction with other species is directly proportional to the product of the masses of the reactants, each raised to a power equal to the coefficient used in the chemical equation. Initially proposed for chemical systems by Cato Guldberg and Peter Waage in 1864 [71] based on previous work by Claude Louis Berthollet [72] and further clarified by Jacobus van’t Hoff in 1877 [73], it also provides a simplified yet effective framework for quantifying biological interactions.

The set of differential equations that corresponds to this interactions is (step-by-step derivation is shown in the Supplementary Materials)

$\begin{array}{cc}\hfill {\displaystyle \frac{d\left[C\right]}{dt}}& =-k_1\left[C\right]\left[Imp\right]+k_{-1}\left[C\right|Imp]+k_{exp}\left[N\right|Exp]\hfill \\ \hfill {\displaystyle \frac{d\left[Imp\right]}{dt}}& =-k_1\left[C\right]\left[Imp\right]+k_{-1}\left[C\right|Imp]+k_{imp}\left[C\right|Imp]\hfill \\ \hfill {\displaystyle \frac{d\left[C\right|Imp]}{dt}}& =k_1\left[C\right]\left[Imp\right]-k_{-1}\left[C\right|Imp]-k_{imp}\left[C\right|Imp]\hfill \\ \hfill {\displaystyle \frac{d\left[N\right]}{dt}}& =-k_2\left[N\right]\left[Exp\right]+k_{-2}\left[N\right|Exp]+k_{imp}\left[C\right|Imp]\hfill \\ \hfill {\displaystyle \frac{d\left[Exp\right]}{dt}}& =-k_2\left[N\right]\left[Exp\right]+k_{-2}\left[N\right|Exp]+k_{exp}\left[N\right|Exp]\hfill \\ \hfill {\displaystyle \frac{d\left[N\right|Exp]}{dt}}& =k_2\left[N\right]\left[Exp\right]-k_{-2}\left[N\right|Exp]-k_{exp}\left[N\right|Exp]\hfill \end{array}$

The numerical solution of this set of ODEs, starting from a given initial species concentration and parameter values (e.g., rate constants for each reaction), provides the temporal evolution of the concentration of all involved species. One important point to take into account is the units of the rate constants: if the time is measured in seconds (s) and the concentrations in mol/L (molar, M), each equation provides the speed of change of the reactants in units of $M/s$. Therefore, for the equations to be consistent in terms of dimensions, the units of the rate constants $k_{-1}$, $k_{-2}$, $k_{imp}$, and $k_{exp}$ have to be in s⁻¹ (first-order reactions), whereas $k_2$ and $k_1$ are measured in (M·s)⁻¹ (second-order reactions in a thermodynamical sense).

Figure 3A–C illustrate three different numerical solutions of the system for an equal set of initial conditions but using the parameter values listed as sets 1, 2, and 3 in Table 1, respectively. In all cases, the concentrations of the species change during the first seconds until they reach a steady state. Both $\left[C\right]$ and $\left[N\right]$ decrease while $\left[C\right|Imp]$ and $\left[N\right|Exp]$ increase (complexes are forming). The enzymes $\left[Imp\right]$ and $\left[Exp\right]$ also decrease due to complex formation. Because of the different parameter values used, three different steady states are reached: In Figure 3A, the final concentration of $\left[C\right]$ in the equilibrium is larger than the concentration of $\left[N\right]$, meaning that most protein is located at the cytoplasm for this value of the rate constant. The opposite scenario is shown in Figure 3C, where the nuclear concentration is now higher. Figure 3B corresponds to a intermediate situation, where the parameters were tuned to obtain an equilibrium condition where the concentrations at the steady state in the nucleus and in the cytoplasm are equivalent. If we imagine a hypothetical situation where our molecule can be seen using immunofluorescence or by direct fusion with a fluorescence protein, case (A) would correspond to a cell that is brighter in the cytoplasm, case (C) would correspond to a cell with a brighter nuclei, and case (B) would correspond to a cell with homogeneous brightness. Note that the different cases were obtained by varying only the values of $k_{imp}$ and $k_{exp}$, meaning that the efficiency of the nuclear import and export processes ultimately controls the balance between the nuclear and cytoplasmic configuration.

6.2. Simplification Using the Mass Conservation Law

The next step in any process of modeling a biological system is to try to simplify and reduce the set of equations. As such, it is often very useful to use the mass conservation law [74], stated by Antoine Lavoisier in 1789 [75].

In the context of biological models of signal transduction, it is quite common to obtain redundant equations where the dynamics of one variable can be simply related to the dynamics of another variable in the system. In this condition, the simplification of the system is direct and results in a more transparent set of equations that is easier to interpret. In this particular case or protein translocation, the mass conservation allowed us to reduce the system of six differential equations to four (step-by-step general application is illustrated in the Supplementary Materials).

As shown in the Supplementary Materials, we obtained three conservation laws for our system. The first conservation corresponds to the fact that the total amount of Smad in all possible configurations (nuclear, cytoplasmic, alone, or in complex) has to remain constant. Two more conservation laws illustrate the fact that the total amount of $Imp$ and $Exp$ proteins has to remain constant. One additional simplification can be introduced, taking into account that at the initial time point, there are no complexes, since we start with all proteins as monomers and monitor their interaction.

In this case, since we often focus on the amount of R-Smads and not in the amount of $Imp$ and $Exp$ proteins, we can use the mass conservation to organize the relationships in such a way that we eliminate the variables $\left[Imp\right]$ and $\left[Exp\right]$, giving us the set of four ODEs:

$\begin{array}{cc}\hfill {\displaystyle \frac{d\left[C\right]}{dt}}& =-k_1\left[C\right](\left[Imp\right]\left(0\right)-\left[C\right|Imp])+\dots \hfill \\ \hfill & +k_{-1}\left[C\right|Imp]+k_{exp}\left[N\right|Exp]\hfill \\ \hfill {\displaystyle \frac{d\left[C\right|Imp]}{dt}}& =k_1\left[C\right](\left[Imp\right]\left(0\right)-\left[C\right|Imp])-\dots \hfill \\ \hfill & -k_{-1}\left[C\right|Imp]-k_{imp}\left[C\right|Imp]\hfill \\ \hfill {\displaystyle \frac{d\left[N\right]}{dt}}& =-k_2\left[N\right](\left[Exp\right]\left(0\right)-\left[N\right|Exp])+\dots \hfill \\ \hfill & +k_{-2}\left[N\right|Exp]+k_{imp}\left[C\right|Imp]\hfill \\ \hfill {\displaystyle \frac{d\left[N\right|Exp]}{dt}}& =k_2\left[N\right](\left[Exp\right]\left(0\right)-\left[N\right|Exp])-\dots \hfill \\ \hfill & -k_{-2}\left[N\right|Exp]-k_{exp}\left[N\right|Exp]\hfill \end{array}$

Numerical integration of this set gives us the temporal evolution of $\left[C\right]$, $\left[N\right]$ and the complexes $\left[C\right|Imp]$ and $\left[N\right|Exp]$. If needed, the temporal evolution of the enzymes can be simply calculated as $\left[Imp\right]\left(t\right)=\left[Imp\right]\left(0\right)-\left[C\right|Imp\left]\right(t)$ and $\left[Exp\right]\left(t\right)=\left[Exp\right]\left(0\right)-\left[N\right|Exp\left]\right(t)$. This means that we can calculate the concentration of all six species over time by solving just four ODEs instead of the six that we obtained initially.

6.3. Nondimensionalization of the Model

Another common approach to increase the predictive capability of a model is to remove the dimensions of the concentrations and the time, expressing both of them as a factor on an appropriate scale for the system under study. In this type of dimensionless system, it is easier to compare parameter values and to interpret the results from the ODEs, since they show relative changes from a maximum to a minimum value instead of absolute values.

In order to make the concentrations dimensionless, we rewrite the model in terms of new variables with no units. This is achieved by dividing each variable by a scaling factor (with concentration units) in such a way that now the concentration of the species $\left[X\right]$ (dimensional) is now represented by $\left[\tilde{X}\right]$ (dimensionless), related via the scaling factor $X_0$ [69]:

(1)$\left[X\right]\left(t\right)=\left[\tilde{X}\right]\left(t\right)·X_0$

The scaling factor for each species is chosen as its maximum value. This way, the dimensionless species concentration $\left[\tilde{X}\right]$ associated falls in the range $[0,1]$. The maximum amount of molecules in the cytoplasm that we can have is $\left[C\right]\left(0\right)$, and, as one molecule of C can only produce one molecule of N, the maximum value for the latter is $\left[C\right]\left(0\right)$ as well (supposing that $\left[N\right]\left(0\right)=0$). In the case of the enzymes, the maximum concentration for each one is also its initial concentration. For the scaling factor of the complexes $C|Imp$ and $N|Exp$, we cannot use their initial concentration as it is zero. Therefore, we must find which is the maximum concentration that they can reach over time (${\left[C\right|Imp]}_{max}$ and ${\left[N\right|Exp]}_{max}$), which happens at a particular time point $t=t_{max}$ (this time value can be different for $\left[C\right|Imp]$ and $\left[N\right|Exp]$):

$\begin{array}{cc}\hfill {\left[C\right|Imp]}_{max}& ={\displaystyle \frac{k_1\left[C\right]\left(t_{max}\right)\left[Imp\right]\left(0\right)}{k_1·\left[C\right]\left(t_{max}\right)+k_{-1}+k_{imp}}}\\ \hfill {\left[N\right|Exp]}_{max}& ={\displaystyle \frac{k_2\left[N\right]\left(t_{max}\right)\left[Exp\right]\left(0\right)}{k_2·\left[N\right]\left(t_{max}\right)+k_{-2}+k_{exp}}}\end{array}$

Next, these expressions can be further simplified by making some assumptions. As the dynamics of the complexes are very fast, $k_2$ and $k_1$ are much larger than the other rate constants. Considering that, we can assume that $k_{-1}=k_{imp}=k_{-2}=k_{exp}=0$; therefore,

$\begin{array}{cc}\hfill {\left[C\right|Imp]}_{max}& \approx \left[Imp\right]\left(0\right)\\ \hfill {\left[N\right|Exp]}_{max}& \approx \left[Exp\right]\left(0\right)\end{array}$

Then, the new dimensionless variables are defined as

$\begin{array}{cc}\hfill \left[\tilde{C}\right]& ={\displaystyle \frac{\left[C\right]}{\left[C\right]\left(0\right)}}\hfill \\ \hfill \left[\widetilde{Imp}\right]& ={\displaystyle \frac{\left[Imp\right]}{\left[Imp\right]\left(0\right)}}\hfill \\ \hfill \left[\widetilde{C|Imp}\right]& ={\displaystyle \frac{\left[C\right|Imp]}{\left[Imp\right]\left(0\right)}}\hfill \\ \hfill \left[\tilde{N}\right]& ={\displaystyle \frac{\left[N\right]}{\left[N\right]\left(0\right)}}\hfill \\ \hfill \left[\widetilde{Exp}\right]& ={\displaystyle \frac{\left[Exp\right]}{\left[Exp\right]\left(0\right)}}\hfill \\ \hfill \left[\widetilde{N|Exp}\right]& ={\displaystyle \frac{\left[N\right|Exp]}{\left[Exp\right]\left(0\right)}}\hfill \end{array}$

Substituting these values into the ODE system, we obtain the new partially dimensionless ODE system:

$\begin{array}{cc}\hfill {\displaystyle \frac{d\left[\tilde{C}\right]}{dt}}& =k_1\left[\tilde{C}\right]\left[Imp\right]\left(0\right)(\left[\widetilde{C|Imp}\right]-1)+\dots \hfill \\ \hfill & +k_{-1}{\displaystyle \frac{\left[Imp\right]\left(0\right)}{\left[C\right]\left(0\right)}}\left[\widetilde{C|Imp}\right]+k_{exp}{\displaystyle \frac{\left[Exp\right]\left(0\right)}{\left[C\right]\left(0\right)}}\left[\widetilde{N|Exp}\right]\hfill \\ \hfill {\displaystyle \frac{d\left[\widetilde{C|Imp}\right]}{dt}}& =k_1\left[\tilde{C}\right]\left[C\right]\left(0\right)(1-\left[\widetilde{C|Imp}\right])-\dots \hfill \\ \hfill & -(k_{-1}+k_{exp})\left[\widetilde{C|Imp}\right]\hfill \\ \hfill {\displaystyle \frac{d\left[\tilde{N}\right]}{dt}}& =k_2\left[\tilde{N}\right]\left[Exp\right]\left(0\right)(\left[\widetilde{N|Exp}\right]-1)+\dots \hfill \\ \hfill & +k_{-2}{\displaystyle \frac{\left[Exp\right]\left(0\right)}{\left[N\right]\left(0\right)}}\left[\widetilde{N|Exp}\right]+k_{imp}{\displaystyle \frac{\left[Imp\right]\left(0\right)}{\left[N\right]\left(0\right)}}\left[\widetilde{C|Imp}\right]\hfill \\ \hfill {\displaystyle \frac{d\left[\widetilde{N|Exp}\right]}{dt}}& =k_2\left[\tilde{N}\right]\left[N\right]\left(0\right)(1-\left[\widetilde{N|Exp}\right])-\dots \hfill \\ \hfill & -(k_{-2}+k_{imp})\left[\widetilde{N|Exp}\right]\hfill \\ \hfill \widetilde{Imp}]\left(t\right)& =1-\left[\widetilde{C|Imp}\right]\left(t\right)\hfill \\ \hfill \widetilde{Exp}]\left(t\right)& =1-\left[\widetilde{N|Exp}\right]\left(t\right)\hfill \end{array}$

It is important to highlight that the mass conservation law in this case is different from the dimensional system. The calculations for this case are not shown, but they are analogous to those in the previous case.

To remove the dimension doe time, a scaling factor $t_0$ is required in order to define a new dimensionless time variable $\tau$:

$\tau ={\displaystyle \frac{t}{t_0}}$

In any system, we have several reaction rate constants that have different values. To obtain the scaling factor, one option is to choose the fastest one, as the dimensionless time is then a factor indicating “how many times the fastest reaction has taken place”. In any case, we can rewrite our ODE system, taking into account that

${\displaystyle \frac{d}{dt}}={\displaystyle \frac{1}{t_0}}{\displaystyle \frac{d}{d\tau }}$

and obtain

$\begin{array}{cc}\hfill {\displaystyle \frac{1}{t_0}}{\displaystyle \frac{d\left[\tilde{C}\right]}{d\tau }}& =k_1\left[\tilde{C}\right]\left[Imp\right]\left(0\right)(\left[\widetilde{C|Imp}\right]-1)+\dots \hfill \\ \hfill & +k_{-1}{\displaystyle \frac{\left[Imp\right]\left(0\right)}{\left[C\right]\left(0\right)}}\left[\widetilde{C|Imp}\right]+k_{exp}{\displaystyle \frac{\left[Exp\right]\left(0\right)}{\left[C\right]\left(0\right)}}\left[\widetilde{N|Exp}\right]\hfill \\ \hfill {\displaystyle \frac{1}{t_0}}{\displaystyle \frac{d\left[\widetilde{C|Imp}\right]}{d\tau }}& =k_1\left[\tilde{C}\right]\left[C\right]\left(0\right)(1-\left[\widetilde{C|Imp}\right])-\dots \hfill \\ \hfill & -(k_{-1}+k_{exp})\left[\widetilde{C|Imp}\right]\hfill \\ \hfill {\displaystyle \frac{1}{t_0}}{\displaystyle \frac{d\left[\tilde{N}\right]}{d\tau }}& =k_2\left[\tilde{N}\right]\left[Exp\right]\left(0\right)(\left[\widetilde{N|Exp}\right]-1)+\dots \hfill \\ \hfill & +k_{-2}{\displaystyle \frac{\left[Exp\right]\left(0\right)}{\left[N\right]\left(0\right)}}\left[\widetilde{N|Exp}\right]+k_{imp}{\displaystyle \frac{\left[Imp\right]\left(0\right)}{\left[N\right]\left(0\right)}}\left[\widetilde{C|Imp}\right]\hfill \\ \hfill {\displaystyle \frac{1}{t_0}}{\displaystyle \frac{d\left[\widetilde{N|Exp}\right]}{d\tau }}& =k_2\left[\tilde{N}\right]\left[N\right]\left(0\right)(1-\left[\widetilde{N|Exp}\right])-\dots \hfill \\ \hfill & -(k_{-2}+k_{imp})\left[\widetilde{N|Exp}\right]\hfill \\ \hfill \widetilde{Imp}]\left(\tau \right)& =1-\left[\widetilde{C|Imp}\right]\left(\tau \right)\hfill \\ \hfill \widetilde{Exp}]\left(\tau \right)& =1-\left[\widetilde{N|Exp}\right]\left(\tau \right)\hfill \end{array}$

In our example, let us suppose that $k_2$ is the fastest rate (which is a biologically plausible assumption). As $k_2$ is measured in (M·s)⁻¹, it cannot be the scaling factor itself because $\tau$ will have dimensions of M·s². Looking at the simplified ODE system using the mass conservation law, we can see that this rate $k_2$ is multiplied by a constant, $\left[Exp\right]\left(0\right)$, measured in M. The product $\left[Exp\right]\left(0\right)·k_2$ is measured in s⁻¹, so we can choose our scaling factor:

$t_0={\displaystyle \frac{1}{k_2·\left[Exp\right]\left(0\right)}}\to \tau =k_2·\left[Exp\right]\left(0\right)·t$

Numerical integration of this dimensionless version for the same parameters is now shown in Figure 4, where the vertical axis now moves from 0 to 1, and the horizontal axis represents multiples of the $t_0$ value.

6.4. Pseudo-Steady-State Approximation

So far, we have an ODE system with 4 coupled equations. Although, in this, case the numerical solution of this system is relatively easy to achieve, this might not be the case for any other system. To further simplify our mathematical model, a common approach is to take advantage of the so-called pseudo-steady-state approximation.

A biologically acceptable approximation is considering that the dynamics of the intermediate complexes are very fast compared with the rest of the species. This means that the translocation processes ($C|Imp$ and $N|Exp$ moving from one compartment to another via diffusion or another active process, plus their disassembly) are slow compared to the association and dissociation reactions ($C|Imp$ and $N|Exp$ formation and disassembly). Mathematically, this translates into $k_{imp}$ and $k_{exp}$ being much smaller than the rest of the rate constants, and thus we can consider that

${\displaystyle \frac{d\left[C\right|Imp]}{dt}}={\displaystyle \frac{d\left[N\right|Exp]}{dt}}=0$

This assumption is known as the Michaelis–Menten approximation [69], in which the concentration of the enzymes and complexes are considered to be in equilibrium while $\left[C\right]$ and $\left[N\right]$ change rapidly. With this mathematical condition, we can then obtain the concentration of the complexes over time as

$\begin{array}{cc}\hfill \left[C\right|Imp\left]\right(t)& ={\displaystyle \frac{k_1\left[C\right]\left(t\right)\left[Imp\right]\left(0\right)}{k_1\left[C\right]\left(t\right)+k_{-1}+k_{imp}}}\hfill \\ \hfill \left[N\right|Exp]\left(t\right)& ={\displaystyle \frac{k_2\left[N\right]\left(t\right)\left[Exp\right]\left(0\right)}{k_2\left[N\right]\left(t\right)+k_{-2}+k_{exp}}}\hfill \end{array}$

With this definition, the new simplified system that gives the concentration of all the species over time is

$\begin{array}{cc}\hfill {\displaystyle \frac{d\left[C\right]}{dt}}& =-k_1\left[C\right]\left[Imp\right]\left(0\right)+\dots \hfill \\ \hfill & +(k_1\left[C\right]+k_{-1}){\displaystyle \frac{k_1\left[C\right]\left[Imp\right]\left(0\right)}{k_1\left[C\right]+k_{-1}+k_{imp}}}+\dots \hfill \\ \hfill & +k_{exp}{\displaystyle \frac{k_2\left[N\right]\left[Exp\right]\left(0\right)}{k_2\left[N\right]+k_{-2}+k_{exp}}}\hfill \\ \hfill {\displaystyle \frac{d\left[N\right]}{dt}}& =-k_2\left[N\right]\left[Exp\right]\left(0\right)+\dots \hfill \\ \hfill & +(k_2\left[N\right]+k_{-2}){\displaystyle \frac{k_2\left[N\right]\left[Exp\right]\left(0\right)}{k_2\left[N\right]+k_{-2}+k_{exp}}}+\dots \hfill \\ \hfill & +k_{imp}{\displaystyle \frac{k_1\left[C\right]\left[Imp\right]\left(0\right)}{k_1\left[C\right]+k_{-1}+k_{imp}}}\hfill \end{array}$

which has just two coupled ODEs. Lastly, using the definition of the Michaelis–Menten constant for each of the two processes (C|$Imp$ and N|$Exp$), we can rewrite the ODEs as

$\begin{array}{cc}\hfill {\displaystyle \frac{d{\left[C\right]}_{MM}}{dt}}& =k_{exp}{\displaystyle \frac{\left[Exp\right]\left(0\right)\left[N\right]}{K_N+\left[N\right]}}-k_{imp}{\displaystyle \frac{\left[Imp\right]\left(0\right)\left[C\right]}{K_C+\left[C\right]}}\\ \hfill {\displaystyle \frac{d{\left[N\right]}_{MM}}{dt}}& =k_{imp}{\displaystyle \frac{\left[Imp\right]\left(0\right)\left[C\right]}{K_C+\left[C\right]}}-k_{exp}{\displaystyle \frac{\left[Exp\right]\left(0\right)\left[N\right]}{K_N+\left[N\right]}}\end{array}$

with:

$\begin{array}{cc}\hfill K_C& ={\displaystyle \frac{k_{imp}+k_{-1}}{k_1}}\\ \hfill K_N& ={\displaystyle \frac{k_{exp}+k_{-2}}{k_2}}\end{array}$

Now, we can compare the initial system with the pseudo-steady-state approximation system that we just obtained. As shown in Figure 5, the approximation is not valid for the values of the rates of set 1 (Table 1), but it is for the values of set 4. The rate constants for the system must fulfill the approximation assumptions in order to be able to use the approximation.

7. Conclusions

Here, we presented a comprehensive revision on how the different compartments in a cell are often utilized as physical barriers to modulate the biochemical processes that shape the response of cells to external stimuli. TGF-$\beta$ is often view as an almost linear cascade, where events at the receptor level correspond to the direct phosphorylation and activation of transcription factors. This oversimplified view leaves out important processes that can fully dictate the response, such as regulation via compartmentalization of the R-Smads, which are actively maintained in and out of the nucleus depending on their activation state. In principle, inactive unphosphorylated R-Smads do not form complexes and do not bind to DNA, so the fact that there is an extra layer that keeps this inactive proteins from entering the nucleus seems unnecessary. But, we know from experience that everything that happens inside the cells often serves many purposes. Therefore, biophysical regulation at least adds an extra safety layer that prevents such important proteins from wrongly activating the expression of their targets when not needed. We believe that this is just one of the many ways that the physical aspects of the cells and organelles influence biological responses, and some are still perhaps waiting to be discovered and studied. For instance, living organisms are open systems by definition, and this property is key to understanding their unique properties, such as adaptation and morphogenesis [76,77,78,79]. Physics and biophysics have already helped to understand many of the most important aspects of how living systems operate, and they will play a central role in future discoveries.

Author Contributions

A.N.F., D.M.-D. and D.G.M. performed this research and wrote the manuscript. D.M.-D. performed numerical simulations. D.G.M. and A.N.F. provided funding. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments

We acknowledge all CBM Severo Ochoa (Centre for Molecular Biology Severo Ochoa, CSIC–UAM) facilities.

Conflicts of Interest

The authors declare no competing interests.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Signaling through TGF-[Forumla omitted. See PDF.] superfamily ligands (adapted from Ref. [29]). In brief, the ligands promote the formation of the receptor complex, which internalizes and phosphorylates R-Smads. After phosphorylation, the R-Smads can form complexes with Smad4 and enter the nuclei, where they bind DNA and activate transcription. Finally, dephosphorylation triggers the disassembly of the complex and the monomers exit the nuclei toward the cytoplasm.

Figure 2. Schematic representation of the nuclear import–export process of a molecule inside a cell. The molecule in the cytoplasm is represented by C and in the nucleus by N. To be imported inside the nucleus, one C molecule must bind one molecule of the import enzyme Imp. Then, the complex [Forumla omitted. See PDF.] can enter the nucleus, where it disassembles into N (former C) and [Forumla omitted. See PDF.]. The export process is analogous, involving the formation of an [Forumla omitted. See PDF.] complex (one molecule of N and one export enzyme [Forumla omitted. See PDF.]) and then the disassembly of this complex in the cytoplasm yielding one molecule of C and one molecule of [Forumla omitted. See PDF.].

Figure 3. Simulation of the ODE system. Three different sets (1, 2, and 3 from Table 1) of rate constants are used for the different cases (A), (B) and (C), respectively. The initial conditions for the species are [Forumla omitted. See PDF.] = [Forumla omitted. See PDF.] = [Forumla omitted. See PDF.] = 0.002 M, [Forumla omitted. See PDF.] = 0.001 M, and [Forumla omitted. See PDF.] = [Forumla omitted. See PDF.] = 0 M in all cases. The concentrations change during the first seconds until they reach a steady state: [Forumla omitted. See PDF.] in equilibrium is bigger (A), equivalent (B), or smaller than [Forumla omitted. See PDF.] in equilibrium. In general, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] tend to decrease initially as [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] form.

Figure 4. Simulation of the ODE system after being made dimensionless. The same rate constants and initial conditions as in Figure 3 were used for all the panels. Note that the concentrations range within [0, 1], and, during the simulation, the curves indicate the present fraction of the concentration of each species in the system (normalized by its scaling factor). In the case of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], the curves show the present fraction of their initial condition. As in Figure 4, the concentrations change during the first seconds until they reach a steady state: [Forumla omitted. See PDF.] in equilibrium is bigger (A), equivalent (B), or smaller than [Forumla omitted. See PDF.] in equilibrium; which for the dimensionless variables (under the conditions for this simulation) means that [Forumla omitted. See PDF.] in equilibrium is equivalent (A), double (B) or greater (C) than [Forumla omitted. See PDF.].

Figure 5. Comparison of the raw model and the pseudo-steady-state approximation model. The concentrations of the C and N molecules are shown as dashed lines for the original ODEs ([Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]) and as solid lines for the pseudo-steady-state approximation (Michaelis–Menten, [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]). In both insets, the initial concentration values are the same as in Figure 3. The rate constants used in (A,B) are sets 1 and 4 (from Table 1), respectively. For case (A), there is a considerable discrepancy between the approximation (solid curves) and the original system (dashed curves). When the rate constants are chosen in such a way that the approximation assumptions are fulfilled (set 4), the steady-state approximation becomes valid and can be used to simplify the analysis of the system.

Table 1

Values of the constant rates for the different sets used in the simulations.

Rate	Set 1	Set 2	Set 3	Set 4
$k_1$	500 (M·s)−1	500 (M·s)−1	500 (M·s)−1	50 (M·s)−1
$k_{-1}$	0.15 s−1	0.15 s−1	0.15 s−1	1.5 s−1
$k_2$	0.03 (M·s)−1	0.036 (M·s)−1	0.04 (M·s)−1	0.03 (M·s)−1
$k_{-2}$	1000 s−1	1000 s−1	1000 s−1	1 s−1
$k_{imp}$	0.1 s−1	0.1 s−1	0.1 s−1	0.3 s−1
$k_{exp}$	0.026 s−1	0.023 s−1	0.01 s−1	0.3 s−1

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/biophysica4040042/s1.

References

1. Heldin, C.H.; Moustakas, A. Signaling Receptors for TGF-beta Family Members. Cold Spring Harb. Perspect. Biol.; 2016; 8, a022053. [DOI: https://dx.doi.org/10.1101/cshperspect.a022053] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27481709]

2. Zhong, Q.; Xiao, X.; Qiu, Y.; Xu, Z.; Chen, C.; Chong, B.; Zhao, X.; Hai, S.; Li, S.; An, Z. et al. Protein posttranslational modifications in health and diseases: Functions, regulatory mechanisms, and therapeutic implications. MedComm; 2023; 4, e261. [DOI: https://dx.doi.org/10.1002/mco2.261] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37143582]

3. Matschinsky, F.M. Editorial: Ebb and Tide of Glucokinase. Int. J. Exp. Diabetes Res.; 2001; 2, pp. 168-171. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12369704]

4. Jetton, T.L.; Shiota, M.; Knobel, S.M.; Piston, D.W.; Cherrington, A.D.; Magnuson, M.A. Substrate-induced Nuclear Export and Peripheral Compartmentalization of Hepatic Glucokinase Correlates with Glycogen Deposition. J. Diabetes Res.; 2001; 2, pp. 173-186. [DOI: https://dx.doi.org/10.1155/EDR.2001.173]

5. Erb, T.J.; Zarzycki, J. A short history of RubisCO: The rise and fall (?) of Nature’s predominant CO₂ fixing enzyme. Curr. Opin. Biotechnol.; 2018; 49, pp. 100-107. [DOI: https://dx.doi.org/10.1016/j.copbio.2017.07.017]

6. Amritkar, K.; Cuevas-Zuviría, B.; Kacar, B. Ancestral structure prediction reveals the conformational impact of the RuBisCO small subunit across time. bioRxiv; 2024; [DOI: https://dx.doi.org/10.1101/2024.06.06.597628]

7. George, C.; Brittain, S.G. A rapid method for quantifying cytoplasmic versus nuclear localization in endogenous peripheral blood leukocytes by conventional flow cytometry. Cytom. Part A; 2017; 91, pp. 351-363. [DOI: https://dx.doi.org/10.1002/cyto.a.23103]

8. Oliver, J.D. ERp57 Functions as a Subunit of Specific Complexes Formed with the ER Lectins Calreticulin and Calnexin. Mol. Biol. Cell; 1999; 10, pp. 2573-2582. [DOI: https://dx.doi.org/10.1091/mbc.10.8.2573]

9. Kimura, T. ERp57 binds competitively to protein disulfide isomerase and calreticulin. Biochem. Biophys. Res. Commun.; 2005; 331, pp. 224-230. [DOI: https://dx.doi.org/10.1016/j.bbrc.2005.03.147]

10. Grillo, C. DNA-binding Activity of the ERp57 C-terminal Domain Is Related to a Redox-dependent Conformational Change. J. Biol. Chem.; 2007; 282, pp. 10299-10310. [DOI: https://dx.doi.org/10.1074/jbc.M700966200]

11. Jeong Sook Kim-Han, K.L.O. Cell Stress Induced by the Parkinsonian Mimetic, 6-Hydroxydopamine, is Concurrent with Oxidation of the Chaperone, ERp57, and Aggresome Formation. Antioxid. Redox Signal.; 2007; 9, pp. 2255-2264. [DOI: https://dx.doi.org/10.1089/ars.2007.1791] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17848102]

12. Grindel, B.J. Tumor necrosis factor-α treatment of HepG2 cells mobilizes a cytoplasmic pool of ERp57/1,25D₃-MARRS to the nucleus. J. Cell. Biochem.; 2011; 112, pp. 2606-2615. [DOI: https://dx.doi.org/10.1002/jcb.23187] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21598303]

13. Scott, M.S. Refining Protein Subcellular Localization. PLoS Comput. Biol.; 2005; 1, e66. [DOI: https://dx.doi.org/10.1371/journal.pcbi.0010066] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16322766]

14. Rajesh Nair, B.R. Protein Subcellular Localization Prediction Using Artificial Intelligence Technology. Methods Mol. Biol.; 2008; 484, pp. 435-463. [DOI: https://dx.doi.org/10.1007/978-1-59745-398-1_27]

15. Mien-Chie Hung, W.L. Protein localization in disease and therapy. J. Cell Sci.; 2011; 124, pp. 3381-3392. [DOI: https://dx.doi.org/10.1242/jcs.089110]

16. Binder, J.X. COMPARTMENTS: Unification and visualization of protein subcellular localization evidence. Database; 2014; 2014, bau012. [DOI: https://dx.doi.org/10.1093/database/bau012]

17. Ruiz-Herrero, T.; Estrada, J.; Guantes, R.; Miguez, D.G. A Tunable Coarse-Grained Model for Ligand-Receptor Interaction. PLoS Comput. Biol.; 2013; 9, e1003274. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1003274]

18. Doldán-Martelli, V.; Míguez, D.G. Synergistic Interaction between Selective Drugs in Cell Populations Models. PLoS ONE; 2015; 10, e0117558. [DOI: https://dx.doi.org/10.1371/journal.pone.0117558]

19. Doldán-Martelli, V.; Guantes, R.; Míguez, D.G. A mathematical model for the rational design of chimeric ligands in selective drug therapies. CPT Pharmacomet. Syst. Pharmacol.; 2013; 2, e26. [DOI: https://dx.doi.org/10.1038/psp.2013.2]

20. Míguez, D.G. The role of asymmetric binding in ligand-receptor systems with 1:2 interaction ratio. Biophys. Chem.; 2010; 148, pp. 74-81. [DOI: https://dx.doi.org/10.1016/j.bpc.2010.02.012]

21. Xu, X.; Zheng, L.; Yuan, Q.; Zhen, G.; Crane, J.L.; Zhou, X.; Cao, X. Transforming growth factor-β in stem cells and tissue homeostasis. Bone Res.; 2018; 6, 2. [DOI: https://dx.doi.org/10.1038/s41413-017-0005-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29423331]

22. Papanikolaou, S.; Vourda, A.; Syggelos, S.; Gyftopoulos, K. Cell Plasticity and Prostate Cancer: The Role of Epithelial–Mesenchymal Transition in Tumor Progression, Invasion, Metastasis and Cancer Therapy Resistance. Cancers; 2021; 13, 2795. [DOI: https://dx.doi.org/10.3390/cancers13112795] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34199763]

23. Derynck, R.; Budi, E.H. Specificity, versatility, and control of TGF-beta family signaling. Sci. Signal.; 2019; 12, eaav5183. [DOI: https://dx.doi.org/10.1126/scisignal.aav5183] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30808818]

24. Deng, Z.; Fan, T.; Xiao, C.; Tian, H.; Zheng, Y.; Li, C.; He, J. TGF-beta signaling in health, disease and therapeutics. Signal Transduct. Target. Ther.; 2024; 9, 61. [DOI: https://dx.doi.org/10.1038/s41392-024-01764-w]

25. Goebel, E.J.; Hart, K.N.; McCoy, J.C.; Thompson, T.B. Structural biology of the TGFβ family. Exp. Biol. Med.; 2019; 244, pp. 1530-1546. [DOI: https://dx.doi.org/10.1177/1535370219880894]

26. Tzavlaki, K.; Moustakas, A. TGF-beta Signaling. Biomolecules; 2020; 10, 487. [DOI: https://dx.doi.org/10.3390/biom10030487]

27. Huang, T.; David, L.; Mendoza, V.; Yang, Y.; De, K.; Sun, L.; Fang, X.; Wrana, J.L.; Hinck, A.P. TGF-β signalling is mediated by two autonomously functioning TBRI:TBRII pairs. EMBO J.; 2011; 30, pp. 1263-1276. [DOI: https://dx.doi.org/10.1038/emboj.2011.54]

28. Mitchell, H.; Choudhury, A.; Pagano, R.E.; Leof, E.B. Ligand-dependent and -independent Transforming Growth Factor-β Receptor Recycling Regulated by Clathrin-mediated Endocytosis and Rab11. Mol. Biol. Cell; 2004; 15, pp. 4166-4178. [DOI: https://dx.doi.org/10.1091/mbc.e04-03-0245]

29. Derynck, R.; Zhang, Y.E. Smad-dependent and Smad-independent pathways in TGF-β family signalling. Nature; 2003; 425, pp. 577-584. [DOI: https://dx.doi.org/10.1038/nature02006]

30. Kim, M.S.; Jin, W. TrkB-Induced Inhibition of R-SMAD/SMAD4 Activation is Essential for TGF-beta-Mediated Tumor Suppressor Activity. Cancers; 2020; 12, 1048. [DOI: https://dx.doi.org/10.3390/cancers12041048]

31. Budi, E.H.; Duan, D.; Derynck, R. Transforming Growth Factor-β Receptors and Smads: Regulatory Complexity and Functional Versatility. Trends Cell Biol.; 2017; 27, pp. 658-672. [DOI: https://dx.doi.org/10.1016/j.tcb.2017.04.005] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28552280]

32. Miyazawa, K.; Miyazono, K. Regulation of TGF-beta Family Signaling by Inhibitory Smads. Cold Spring Harb. Perspect. Biol.; 2017; 9, a022095. [DOI: https://dx.doi.org/10.1101/cshperspect.a022095] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27920040]

33. Aykul, S.; Maust, J.; Thamilselvan, V.; Floer, M.; Martinez-Hackert, E. Smad2/3 Activation Regulates Smad1/5/8 Signaling via a Negative Feedback Loop to Inhibit 3T3-L1 Adipogenesis. Int. J. Mol. Sci.; 2021; 22, 8472. [DOI: https://dx.doi.org/10.3390/ijms22168472] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34445177]

34. Baburajendran, N.; Jauch, R.; Tan, C.Y.Z.; Narasimhan, K.; Kolatkar, P.R. Structural basis for the cooperative DNA recognition by Smad4 MH1 dimers. Nucleic Acids Res.; 2011; 39, pp. 8213-8222. [DOI: https://dx.doi.org/10.1093/nar/gkr500]

35. Xin, X.; Cheng, X.; Zeng, F.; Xu, Q.; Hou, L. The Role of TGF-beta/SMAD Signaling in Hepatocellular Carcinoma: From Mechanism to Therapy and Prognosis. Int. J. Biol. Sci.; 2024; 20, pp. 1436-1451. [DOI: https://dx.doi.org/10.7150/ijbs.89568]

36. Hata, A.; Chen, Y.G. TGF-beta Signaling from Receptors to Smads. Cold Spring Harb. Perspect. Biol.; 2016; 8, a022061. [DOI: https://dx.doi.org/10.1101/cshperspect.a022061]

37. Xu, F.; Liu, C.; Zhou, D.; Zhang, L. TGF-beta/SMAD Pathway and Its Regulation in Hepatic Fibrosis. J. Histochem. Cytochem.; 2016; 64, pp. 157-167. [DOI: https://dx.doi.org/10.1369/0022155415627681]

38. Bai, H.; Lester, G.M.; Petishnok, L.; Dean, D. Cytoplasmic transport and nuclear import of plasmid DNA. Biosci. Rep.; 2017; 37, BSR20160616. [DOI: https://dx.doi.org/10.1042/BSR20160616]

39. Schmierer, B.; Hill, C.S. Kinetic Analysis of Smad Nucleocytoplasmic Shuttling Reveals a Mechanism for Transforming Growth Factor β-Dependent Nuclear Accumulation of Smads. Mol. Cell. Biol.; 2005; 25, pp. 9845-9858. [DOI: https://dx.doi.org/10.1128/MCB.25.22.9845-9858.2005]

40. Hill, C.S. Nucleocytoplasmic shuttling of Smad proteins. Cell Res.; 2009; 19, pp. 36-46. [DOI: https://dx.doi.org/10.1038/cr.2008.325]

41. Otsuka, S.; Ellenberg, J. Mechanisms of nuclear pore complex assembly—Two different ways of building one molecular machine. FEBS Lett.; 2018; 592, pp. 475-488. [DOI: https://dx.doi.org/10.1002/1873-3468.12905] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29119545]

42. Paci, G.; Caria, J.; Lemke, E.A. Cargo transport through the nuclear pore complex at a glance. J. Cell Sci.; 2021; 134, jcs247874. [DOI: https://dx.doi.org/10.1242/jcs.247874] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33495357]

43. Timney, B.L.; Raveh, B.; Mironska, R.; Trivedi, J.M.; Kim, S.J.; Russel, D.; Wente, S.R.; Sali, A.; Rout, M.P. Simple rules for passive diffusion through the nuclear pore complex. J. Cell Biol.; 2016; 215, pp. 57-76. [DOI: https://dx.doi.org/10.1083/jcb.201601004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27697925]

44. Oka, M.; Yoneda, Y. Importin α: Functions as a nuclear transport factor and beyond. Proc. Jpn. Acad. Ser. B; 2018; 94, pp. 259-274. [DOI: https://dx.doi.org/10.2183/pjab.94.018] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30078827]

45. Yang, Y.; Guo, L.; Chen, L.; Gong, B.; Jia, D.; Sun, Q. Nuclear transport proteins: Structure, function and disease relevance. Signal Transduct. Target. Ther.; 2023; 8, 425. [DOI: https://dx.doi.org/10.1038/s41392-023-01649-4]

46. Kurisaki, A.; Kose, S.; Yoneda, Y.; Heldin, C.H.; Moustakas, A. Transforming Growth Factor-β Induces Nuclear Import of Smad3 in an Importin-β1 and Ran-dependent Manner. Mol. Biol. Cell; 2001; 12, pp. 1079-1091. [DOI: https://dx.doi.org/10.1091/mbc.12.4.1079]

47. Chen, X.; Xu, L. Mechanism and Regulation of Nucleocytoplasmic Trafficking of Smad. Cell Biosci.; 2011; 1, 40. [DOI: https://dx.doi.org/10.1186/2045-3701-1-40]

48. Kırlı, K.; Karaca, S.; Dehne, H.J.; Samwer, M.; Pan, K.T.; Lenz, C.; Urlaub, H.; Görlich, D. A deep proteomics perspective on CRM1-mediated nuclear export and nucleocytoplasmic partitioning. eLife; 2015; 4, e11466. [DOI: https://dx.doi.org/10.7554/eLife.11466]

49. Abdel Mouti, M.; Pauklin, S. TGFB1/INHBA Homodimer/Nodal-SMAD2/3 Signaling Network: A Pivotal Molecular Target in PDAC Treatment. Mol. Ther.; 2021; 29, pp. 920-936. [DOI: https://dx.doi.org/10.1016/j.ymthe.2021.01.002]

50. Watanabe, M.; Masuyama, N.; Fukuda, M.; Nishida, E. Regulation of intracellular dynamics of Smad4 by its leucine-rich nuclear export signal. EMBO Rep.; 2000; 1, pp. 176-182. [DOI: https://dx.doi.org/10.1093/embo-reports/kvd029]

51. Dai, F.; Lin, X.; Chang, C.; Feng, X.H. Nuclear Export of Smad2 and Smad3 by RanBP3 Facilitates Termination of TGF-β Signaling. Dev. Cell; 2009; 16, pp. 345-357. [DOI: https://dx.doi.org/10.1016/j.devcel.2009.01.022] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19289081]

52. Xu, L.; Kang, Y.; Çöl, S.; Massagué, J. Smad2 Nucleocytoplasmic Shuttling by Nucleoporins CAN/Nup214 and Nup153 Feeds TGFβ Signaling Complexes in the Cytoplasm and Nucleus. Mol. Cell; 2002; 10, pp. 271-282. [DOI: https://dx.doi.org/10.1016/S1097-2765(02)00586-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12191473]

53. Kofler, M.; Kapus, A. Nuclear Import and Export of YAP and TAZ. Cancers; 2023; 15, 4956. [DOI: https://dx.doi.org/10.3390/cancers15204956] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37894323]

54. Kurisaki, A.; Kurisaki, K.; Kowanetz, M.; Sugino, H.; Yoneda, Y.; Heldin, C.H.; Moustakas, A. The Mechanism of Nuclear Export of Smad3 Involves Exportin 4 and Ran. Mol. Cell. Biol.; 2006; 26, pp. 1318-1332. [DOI: https://dx.doi.org/10.1128/MCB.26.4.1318-1332.2006]

55. Song, R.; Tikoo, S.; Jain, R.; Pinello, N.; Au, A.Y.M.; Nagarajah, R.; Porse, B.; Rasko, J.E.J.; Wong, J.J.-L. Dynamic intron retention modulates gene expression in the monocytic differentiation pathway. Immunology; 2022; 165, pp. 274-286. [DOI: https://dx.doi.org/10.1111/imm.13435]

56. Pierreux, C.E.; Nicolás, F.J.; Hill, C.S. Transforming Growth Factor β-Independent Shuttling of Smad4 between the Cytoplasm and Nucleus. Mol. Cell. Biol.; 2000; 20, pp. 9041-9054. [DOI: https://dx.doi.org/10.1128/MCB.20.23.9041-9054.2000]

57. Zou, M.L.; Chen, Z.H.; Teng, Y.Y.; Liu, S.Y.; Jia, Y.; Zhang, K.W.; Sun, Z.L.; Wu, J.J.; Yuan, Z.D.; Feng, Y. et al. The Smad Dependent TGF-beta and BMP Signaling Pathway in Bone Remodeling and Therapies. Front. Mol. Biosci.; 2021; 8, 593310. [DOI: https://dx.doi.org/10.3389/fmolb.2021.593310]

58. Warmflash, A.; Zhang, Q.; Sorre, B.; Vonica, A.; Brivanlou, A.H. Dynamics of TGF-beta signaling reveal adaptive and pulsatile behaviors reflected in the nuclear localization of transcription factor Smad4. Proc. Natl. Acad. Sci. USA; 2012; 109, pp. E1947-E1956. [DOI: https://dx.doi.org/10.1073/pnas.1207607109]

59. Li, Y.; Luo, W.; Yang, W. Nuclear Transport and Accumulation of Smad Proteins Studied by Single-Molecule Microscopy. Biophys. J.; 2018; 114, pp. 2243-2251. [DOI: https://dx.doi.org/10.1016/j.bpj.2018.03.018]

60. Nicolás, F.J.; De Bosscher, K.; Schmierer, B.; Hill, C.S. Analysis of Smad nucleocytoplasmic shuttling in living cells. J. Cell Sci.; 2004; 117, pp. 4113-4125. [DOI: https://dx.doi.org/10.1242/jcs.01289]

61. Zi, Z.; Feng, Z.; Chapnick, D.A.; Dahl, M.; Deng, D.; Klipp, E.; Moustakas, A.; Liu, X. Quantitative analysis of transient and sustained transforming growth factor-β signaling dynamics. Mol. Syst. Biol.; 2011; 7, 492. [DOI: https://dx.doi.org/10.1038/msb.2011.22] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21613981]

62. Sorre, B.; Warmflash, A.; Brivanlou, A.; Siggia, E. Encoding of Temporal Signals by the TGF-β Pathway and Implications for Embryonic Patterning. Dev. Cell; 2014; 30, pp. 334-342. [DOI: https://dx.doi.org/10.1016/j.devcel.2014.05.022] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25065773]

63. Frick, C.L.; Yarka, C.; Nunns, H.; Goentoro, L. Sensing relative signal in the Tgf-β/Smad pathway. Proc. Natl. Acad. Sci. USA; 2017; 114, pp. E2975-E2982. [DOI: https://dx.doi.org/10.1073/pnas.1611428114] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28320972]

64. Kretzschmar, M.; Doody, J.; Timokhina, I.; Massague, J. A mechanism of repression of TGFbeta/Smad signaling by oncogenic Ras. Genes Dev.; 1999; 13, pp. 804-816. [DOI: https://dx.doi.org/10.1101/gad.13.7.804]

65. Míguez, D.G.; Gil-Guiñón, E.; Pons, S.; Martí, E. Smad2 and Smad3 cooperate and antagonize simultaneously in vertebrate neurogenesis. J. Cell Sci.; 2013; 126, pp. 5335-5343. [DOI: https://dx.doi.org/10.1242/jcs.130435]

66. Sanz, A.L.; Míguez, D.G. Dual R-Smads interplay in the regulation of vertebrate neurogenesis. Neurogenesis; 2014; 1, e29529. [DOI: https://dx.doi.org/10.4161/neur.29529]

67. Schmierer, B.; Tournier, A.L.; Bates, P.A.; Hill, C.S. Mathematical modeling identifies Smad nucleocytoplasmic shuttling as a dynamic signal-interpreting system. Proc. Natl. Acad. Sci. USA; 2008; 105, pp. 6608-6613. [DOI: https://dx.doi.org/10.1073/pnas.0710134105]

68. Khatibi, S.; Zhu, H.J.; Wagner, J.; Tan, C.W.; Burgess, A.W. Mathematical model of TGF- beta signalling: Feedback coupling is consistent with signal switching. BMC Syst. Biol.; 2017; 11, 48. [DOI: https://dx.doi.org/10.1186/s12918-017-0421-5]

69. Chen, W.W.; Niepel, M.; Sorger, P.K. Classic and contemporary approaches to modeling biochemical reactions. Genes Dev.; 2010; 24, pp. 1861-1875. [DOI: https://dx.doi.org/10.1101/gad.1945410]

70. Hangos, K.M.; Szederkényi, G. Mass action realizations of reaction kinetic system models on various time scales. J. Phys. Conf. Ser.; 2011; 268, 012009. [DOI: https://dx.doi.org/10.1088/1742-6596/268/1/012009]

71. Guldberg, C.M.; Waage, P. Ueber die chemische Affinität. 1. Einleitung. J. Prakt. Chem.; 1879; 19, pp. 69-114. [DOI: https://dx.doi.org/10.1002/prac.18790190111]

72. Berthollet, C.L. Essai de Statique Chimique; A Paris: Chez Firmin Didot … Nabu Press: New York, NY, USA, 1803.

73. Van’T Hoff, J.H. Die Grenzebene, ein Beitrag zur Kenntniss der Esterbildung. Ber. Dtsch. Chem. Ges.; 1877; 10, pp. 669-678. [DOI: https://dx.doi.org/10.1002/cber.187701001185]

74. Weiss, L. A simple demonstration of the law of conservation of mass. J. Chem. Educ.; 1943; 20, 579. [DOI: https://dx.doi.org/10.1021/ed020p579]

75. Schummer, J. The Philosophy of Chemistry. Philosophy of Chemistry: Synthesis of a New Discipline; Baird, D.; Scerri, E.; McIntyre, L. Springer: Dordrecht, The Netherlands, 2006; pp. 19-39. [DOI: https://dx.doi.org/10.1007/1-4020-3261-7_2]

76. Míguez, D.G.; Satnoianu, R.; Muñuzuri, A. Experimental steady pattern formation in reaction-diffusion-advection systems. Phys. Rev. E; 2006; 73, 25201. [DOI: https://dx.doi.org/10.1103/PhysRevE.73.025201] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16605385]

77. Míguez, D.G.; Izús, G.; Muñuzuri, A. Robustness and stability of flow-and-diffusion structures. Phys. Rev. E; 2006; 73, 016207. [DOI: https://dx.doi.org/10.1103/PhysRevE.73.016207]

78. Pierre Madl, P.R. Quantum Electrodynamics Coherence and Hormesis: Foundations of Quantum Biology. Int. J. Mol. Sci.; 2023; 24, 14003. [DOI: https://dx.doi.org/10.3390/ijms241814003]

79. Renati, P. Reframing Some Methods in Life Sciences: Relationship and Causation in Living Matter. Cut. Edge Res. Biol.; 2023; 8, pp. 127-194. [DOI: https://dx.doi.org/10.9734/bpi/cerb/v8/19273d]