Technologies for whole-cell modeling: Genome-wide reconstruction of a cell in silico.

With advances in high-throughput, large-scale in vivo measurement and genome modification techniques at the single-nucleotide level, there is an increasing demand for the development of new technologies for the flexible design and control of cellular systems. Computer-aided design is a powerful tool to design new cells. Whole-cell modeling aims to integrate various cellular subsystems, determine their interactions and cooperative mechanisms, and predict comprehensive cellular behaviors by computational simulations on a genome-wide scale. It has been applied to prokaryotes, yeasts, and higher eukaryotic cells, and utilized in a wide range of applications, including production of valuable substances, drug discovery, and controlled differentiation. Whole-cell modeling, consisting of several thousand elements with diverse scales and properties, requires innovative model construction, simulation, and analysis techniques. Furthermore, whole-cell modeling has been extended to multiple scales, including high-resolution modeling at the single-nucleotide and single-amino acid levels and multicellular modeling of tissues and organs. This review presents an overview of the current state of whole-cell modeling, discusses the novel computational and experimental technologies driving it, and introduces further developments toward multihierarchical modeling on a whole-genome scale.

eGFRD in all dimensions.

Biochemical reactions often occur at low copy numbers but at once in crowded and diverse environments. Space and stochasticity therefore play an essential role in biochemical networks. Spatial-stochastic simulations have become a prominent tool for understanding how stochasticity at the microscopic level influences the macroscopic behavior of such systems. While particle-based models guarantee the level of detail necessary to accurately describe the microscopic dynamics at very low copy numbers, the algorithms used to simulate them typically imply trade-offs between computational efficiency and biochemical accuracy. eGFRD (enhanced Green's Function Reaction Dynamics) is an exact algorithm that evades such trade-offs by partitioning the N-particle system into M ≤ N analytically tractable one- and two-particle systems; the analytical solutions (Green's functions) then are used to implement an event-driven particle-based scheme that allows particles to make large jumps in time and space while retaining access to their state variables at arbitrary simulation times. Here we present "eGFRD2," a new eGFRD version that implements the principle of eGFRD in all dimensions, thus enabling efficient particle-based simulation of biochemical reaction-diffusion processes in the 3D cytoplasm, on 2D planes representing membranes, and on 1D elongated cylinders representative of, e.g., cytoskeletal tracks or DNA; in 1D, it also incorporates convective motion used to model active transport. We find that, for low particle densities, eGFRD2 is up to 6 orders of magnitude faster than conventional Brownian dynamics. We exemplify the capabilities of eGFRD2 by simulating an idealized model of Pom1 gradient formation, which involves 3D diffusion, active transport on microtubules, and autophosphorylation on the membrane, confirming recent experimental and theoretical results on this system to hold under genuinely stochastic conditions.

The Berg-Purcell limit revisited.

Biological systems often have to measure extremely low concentrations of chemicals with high precision. When dealing with such small numbers of molecules, the inevitable randomness of physical transport processes and binding reactions will limit the precision with which measurements can be made. An important question is what the lower bound on the noise would be in such measurements. Using the theory of diffusion-influenced reactions, we derive an analytical expression for the precision of concentration estimates that are obtained by monitoring the state of a receptor to which a diffusing ligand can bind. The variance in the estimate consists of two terms, one resulting from the intrinsic binding kinetics and the other from the diffusive arrival of ligand at the receptor. The latter term is identical to the fundamental limit derived by Berg and Purcell (Biophys. J., 1977), but disagrees with a more recent expression by Bialek and Setayeshgar. Comparing the theoretical predictions against results from particle-based simulations confirms the accuracy of the resulting expression and reaffirms the fundamental limit established by Berg and Purcell.

Transomics2cytoscape: an automated software for interpretable 2.5-dimensional visualization of trans-omic networks.

Biochemical network visualization is one of the essential technologies for mechanistic interpretation of omics data. In particular, recent advances in multi-omics measurement and analysis require the development of visualization methods that encompass multiple omics data. Visualization in 2.5 dimension (2.5D visualization), which is an isometric view of stacked X-Y planes, is a convenient way to interpret multi-omics/trans-omics data in the context of the conventional layouts of biochemical networks drawn on each of the stacked omics layers. However, 2.5D visualization of trans-omics networks is a state-of-the-art method that primarily relies on time-consuming human efforts involving manual drawing. Here, we present an R Bioconductor package 'transomics2cytoscape' for automated visualization of 2.5D trans-omics networks. We confirmed that transomics2cytoscape could be used for rapid visualization of trans-omics networks presented in published papers within a few minutes. Transomics2cytoscape allows for frequent update/redrawing of trans-omics networks in line with the progress in multi-omics/trans-omics data analysis, thereby enabling network-based interpretation of multi-omics data at each research step. The transomics2cytoscape source code is available at https://github.com/ecell/transomics2cytoscape .

An architecture for collaboration in systems biology at the age of the Metaverse.

As the current state of the Metaverse is largely driven by corporate interests, which may not align with scientific goals and values, academia should play a more active role in its development. Here, we present the challenges and solutions for building a Metaverse that supports systems biology research and collaboration. Our solution consists of two components: Kosmogora, a server ensuring biological data access, traceability, and integrity in the context of a highly collaborative environment such as a metaverse; and ECellDive, a virtual reality application to explore, interact, and build upon the data managed by Kosmogora. We illustrate the synergy between the two components by visualizing a metabolic network and its flux balance analysis. We also argue that the Metaverse of systems biology will foster closer communication and cooperation between experimentalists and modelers in the field.

Fragilities caused by dosage imbalance in regulation of the budding yeast cell cycle.

Cells can maintain their functions despite fluctuations in intracellular parameters, such as protein activities and gene expression levels. This commonly observed biological property of cells is called robustness. On the other hand, these parameters have different limitations, each reflecting the property of the subsystem containing the parameter. The budding yeast cell cycle is quite fragile upon overexpression of CDC14, but is robust upon overexpression of ESP1. The gene products of both CDC14 and ESP1 are regulated by 1ratio1 binding with their inhibitors (Net1 and Pds1), and a mathematical model predicts the extreme fragility of the cell cycle upon overexpression of CDC14 and ESP1 caused by dosage imbalance between these genes. However, it has not been experimentally shown that dosage imbalance causes fragility of the cell cycle. In this study, we measured the quantitative genetic interactions of these genes by performing combinatorial "genetic tug-of-war" experiments. We first showed experimental evidence that dosage imbalance between CDC14 and NET1 causes fragility. We also showed that fragility arising from dosage imbalance between ESP1 and PDS1 is masked by CDH1 and CLB2. The masking function of CLB2 was stabilization of Pds1 by its phosphorylation. We finally modified Chen's model according to our findings. We thus propose that dosage imbalance causes fragility in biological systems.

A comprehensive molecular interaction map of the budding yeast cell cycle.

With the accumulation of data on complex molecular machineries coordinating cell-cycle dynamics, coupled with its central function in disease patho-physiologies, it is becoming increasingly important to collate the disparate knowledge sources into a comprehensive molecular network amenable to systems-level analyses. In this work, we present a comprehensive map of the budding yeast cell-cycle, curating reactions from ∼600 original papers. Toward leveraging the map as a framework to explore the underlying network architecture, we abstract the molecular components into three planes--signaling, cell-cycle core and structural planes. The planar view together with topological analyses facilitates network-centric identification of functions and control mechanisms. Further, we perform a comparative motif analysis to identify around 194 motifs including feed-forward, mutual inhibitory and feedback mechanisms contributing to cell-cycle robustness. We envisage the open access, comprehensive cell-cycle map to open roads toward community-based deeper understanding of cell-cycle dynamics.

Signaling activations through G-protein-coupled-receptor aggregations.

Eukaryotic cells transmit extracellular signal information to cellular interiors through the formation of a ternary complex made up of a ligand (or agonist), G-protein, and G-protein-coupled receptor (GPCR). Previously formalized theories of ternary complex formation have mainly assumed that observable states of receptors can only take the form of monomers. Here, we propose a multiary complex model of GPCR signaling activations via the vector representation of various unobserved aggregated receptor states. Our results from model simulations imply that receptor aggregation processes can govern cooperative effects in a regime inaccessible by previous theories. In particular, we show how the affinity of ligand-receptor binding can be largely varied by various oligomer formations in the low concentration range of G-protein stimulus.

Simulation of live-cell imaging system reveals hidden uncertainties in cooperative binding measurements.

We propose a computational method to quantitatively evaluate the systematic uncertainties that arise from undetectable sources in biological measurements using live-cell imaging techniques. We then demonstrate this method in measuring the biological cooperativity of molecular binding networks, in particular, ligand molecules binding to cell-surface receptor proteins. Our results show how the nonstatistical uncertainties lead to invalid identifications of the measured cooperativity. Through this computational scheme, the biological interpretation can be more objectively evaluated and understood under a specific experimental configuration of interest.

Cooperativity transitions driven by higher-order oligomer formations in ligand-induced receptor dimerization.

While cooperativity in ligand-induced receptor dimerization has been linked with receptor-receptor couplings via minimal representations of physical observables, effects arising from higher-order oligomer, e.g., trimer and tetramer, formations of unobserved receptors have received less attention. Here we propose a dimerization model of ligand-induced receptors in multivalent form representing physical observables under basis vectors of various aggregated receptor states. Our simulations of multivalent models not only reject Wofsy-Goldstein parameter conditions for cooperativity, but show that higher-order oligomer formations can shift cooperativity from positive to negative.

Surface reaction-diffusion kinetics on lattice at the microscopic scale.

Microscopic models of reaction-diffusion processes on the cell membrane can link local spatiotemporal effects to macroscopic self-organized patterns often observed on the membrane. Simulation schemes based on the microscopic lattice method (MLM) can model these processes at the microscopic scale by tracking individual molecules, represented as hard spheres, on fine lattice voxels. Although MLM is simple to implement and is generally less computationally demanding than off-lattice approaches, its accuracy and consistency in modeling surface reactions have not been fully verified. Using the Spatiocyte scheme, we study the accuracy of MLM in diffusion-influenced surface reactions. We derive the lattice-based bimolecular association rates for two-dimensional (2D) surface-surface reaction and one-dimensional (1D) volume-surface adsorption according to the Smoluchowski-Collins-Kimball model and random walk theory. We match the time-dependent rates on lattice with off-lattice counterparts to obtain the correct expressions for MLM parameters in terms of physical constants. The expressions indicate that the voxel size needs to be at least 0.6% larger than the molecule to accurately simulate surface reactions on triangular lattice. On square lattice, the minimum voxel size should be even larger, at 5%. We also demonstrate the ability of MLM-based schemes such as Spatiocyte to simulate a reaction-diffusion model that involves all dimensions: three-dimensional (3D) diffusion in the cytoplasm, 2D diffusion on the cell membrane, and 1D cytoplasm-membrane adsorption. With the model, we examine the contribution of the 2D reaction pathway to the overall reaction rate at different reactant diffusivity, reactivity, and concentrations.

A mathematical model for the Kai-protein-based chemical oscillator and clock gene expression rhythms in cyanobacteria.

In the cyanobacterium, Synechococcus elongatus, most promoters are regulated by a circadian clock under continuous light (LL) conditions. Nevertheless, the basic circadian oscillation is primarily generated by alternating KaiC phosphorylation/dephosphorylation reactions at the posttranslational level. Indeed, the KaiC phosphorylation cycle was recently reconstituted in vitro by incubating KaiA, KaiB, and KaiC proteins with ATP. However, the molecular dynamics of this chemical oscillation and the mechanism that drives the circadian transcription/translation rhythms remain unknown. In this report, the KaiC phosphorylation cycle and the gene regulatory network in the cyanobacterial circadian system have been modeled. The model reproduces the robust KaiC phosphorylation cycle in the absence of de novo gene expression as is observed in vitro, as well as its coupling to transcriptional/translational feedback in LL conditions in vivo. Moreover, the model is consistent with most previous experiments, including various combinations of genetic knockout or overexpression of kai genes. It also predicts that multiple KaiC phosphorylation states and dynamic Kai protein interactions may be required for the cyanobacterial circadian system.

Transient Acceleration of Epidermal Growth Factor Receptor Dynamics Produces Higher-Order Signaling Clusters.

Cell signaling depends on spatiotemporally regulated molecular interactions. Although the movements of signaling proteins have been analyzed with various technologies, how spatial dynamics influence the molecular interactions that transduce signals is unclear. Here, we developed a single-molecule method to analyze the spatiotemporal coupling between motility, clustering, and signaling. The analysis was performed with the epidermal growth factor receptor (EGFR), which triggers signaling through its dimerization and phosphorylation after association with EGF. Our results show that the few EGFRs isolated in membrane subdomains were released by an EGF-dependent increase in their diffusion area, facilitating molecular associations and producing immobile clusters. Using a two-color single-molecule analysis, we found that the EGF-induced state transition alters the properties of the immobile clusters, allowing them to interact for extended periods with the cytoplasmic protein, GRB2. Our study reveals a novel correlation between this molecular interaction and its mesoscale dynamics, providing the initial signaling node.

Density imaging of heterochromatin in live cells using orientation-independent-DIC microscopy.

In eukaryotic cells, highly condensed inactive/silenced chromatin has long been called "heterochromatin." However, recent research suggests that such regions are in fact not fully transcriptionally silent and that there exists only a moderate access barrier to heterochromatin. To further investigate this issue, it is critical to elucidate the physical properties of heterochromatin such as its total density in live cells. Here, using orientation-independent differential interference contrast (OI-DIC) microscopy, which is capable of mapping optical path differences, we investigated the density of the total materials in pericentric foci, a representative heterochromatin model, in live mouse NIH3T3 cells. We demonstrated that the total density of heterochromatin (208 mg/ml) was only 1.53-fold higher than that of the surrounding euchromatic regions (136 mg/ml) while the DNA density of heterochromatin was 5.5- to 7.5-fold higher. We observed similar minor differences in density in typical facultative heterochromatin, the inactive human X chromosomes. This surprisingly small difference may be due to that nonnucleosomal materials (proteins/RNAs) (∼120 mg/ml) are dominant in both chromatin regions. Monte Carlo simulation suggested that nonnucleosomal materials contribute to creating a moderate access barrier to heterochromatin, allowing minimal protein access to functional regions. Our OI-DIC imaging offers new insight into the live cellular environments.

The physical size of transcription factors is key to transcriptional regulation in chromatin domains.

Genetic information, which is stored in the long strand of genomic DNA as chromatin, must be scanned and read out by various transcription factors. First, gene-specific transcription factors, which are relatively small (∼50 kDa), scan the genome and bind regulatory elements. Such factors then recruit general transcription factors, Mediators, RNA polymerases, nucleosome remodellers, and histone modifiers, most of which are large protein complexes of 1-3 MDa in size. Here, we propose a new model for the functional significance of the size of transcription factors (or complexes) for gene regulation of chromatin domains. Recent findings suggest that chromatin consists of irregularly folded nucleosome fibres (10 nm fibres) and forms numerous condensed domains (e.g., topologically associating domains). Although the flexibility and dynamics of chromatin allow repositioning of genes within the condensed domains, the size exclusion effect of the domain may limit accessibility of DNA sequences by transcription factors. We used Monte Carlo computer simulations to determine the physical size limit of transcription factors that can enter condensed chromatin domains. Small gene-specific transcription factors can penetrate into the chromatin domains and search their target sequences, whereas large transcription complexes cannot enter the domain. Due to this property, once a large complex binds its target site via gene-specific factors it can act as a 'buoy' to keep the target region on the surface of the condensed domain and maintain transcriptional competency. This size-dependent specialization of target-scanning and surface-tethering functions could provide novel insight into the mechanisms of various DNA transactions, such as DNA replication and repair/recombination.

Evaluation of the lower protein limit in the budding yeast Saccharomyces cerevisiae using TIPI-gTOW.

BACKGROUND: Identifying permissible limits of intracellular parameters such as protein expression provides important information for examining robustness. In this study, we used the TEV protease-mediated induction of protein instability (TIPI) in combination with the genetic Tug-of-War (gTOW) to develop a method to measure the lower limit of protein level. We first tested the feasibility of this method using ADE2 as a marker and then analyzed some cell cycle regulators to reveal genetic interactions. RESULTS: Using TIPI-gTOW, we successfully constructed a strain in which GFP-(TDegF)Ade2 was expressed at the lower limit, just sufficient to support cellular growth under the -Ade condition by accelerating degradation by TEV protease. We also succeeded in constructing a strain in which the minimal level of GFP-(TDegF)Cdc20 was expressed by TIPI-gTOW. Using this strain, we studied genetic interactions between cell cycle regulators and CDC20, and the result was highly consistent with the previously identified interactions. Comparison of the experimental data with predictions of a mathematical model revealed some interactions that were not implemented into the current model. CONCLUSIONS: TIPI-gTOW is useful for estimating changes in the lower limit of a protein under different conditions, such as different genetic backgrounds and environments. TIPI-gTOW is also useful for analyzing genetic interactions of essential genes whose deletion mutants cannot be obtained.

A quantitative model of ERK MAP kinase phosphorylation in crowded media.

Cytoplasm contains a large number of macromolecules at extremely high densities. How this striking nature of intracellular milieu called macromolecular crowding affects intracellular signaling remains uncharacterized. Here, we examined the effect of macromolecular crowding on ERK MAPK phosphorylation by MEK MAPKK. Addition of polyethylene glycol-6000 (PEG-6000) as a crowder to mimic intracellular environments, elicited a biphasic response to the overall ERK phosphorylation rate. Furthermore, probability of processive phosphorylation (processivity) of tyrosine and threonine residues within the activation loop on ERK increased non-linearly for increasing PEG-6000 concentration. Based on the experimental data, we developed for the first time a mathematical model integrating all of the effects of thermodynamic activity, viscosity, and processivity in crowded media, and found that ERK phosphorylation is transition-state-limited reaction. The mathematical model allows accurate estimation of the effects of macromolecular crowding on a wide range of reaction kinetics, from transition-state-limited to diffusion-limited reactions.

Flexible and dynamic nucleosome fiber in living mammalian cells.

Genomic DNA is organized three dimensionally within cells as chromatin and is searched and read by various proteins by an unknown mechanism; this mediates diverse cell functions. Recently, several pieces of evidence, including our cryomicroscopy and synchrotron X-ray scattering analyses, have demonstrated that chromatin consists of irregularly folded nucleosome fibers without a 30-nm chromatin fiber (i.e., a polymer melt-like structure). This melt-like structure implies a less physically constrained and locally more dynamic state, which may be crucial for protein factors to scan genomic DNA. Using a combined approach of fluorescence correlation spectroscopy, Monte Carlo computer simulations, and single nucleosome imaging, we demonstrated the flexible and dynamic nature of the nucleosome fiber in living mammalian cells. We observed local nucleosome fluctuation (~50 nm movement/30 ms) caused by Brownian motion. Our in vivo/in silico results suggest that local nucleosome dynamics facilitate chromatin accessibility and play a critical role in the scanning of genome information.