Biomolecular functions emerge from the molecular interactions taking place in cellular environments. Understanding each component’s role in driving cell function poses an immense challenge. For a long time, experimental techniques have been our main window into the cellular environment. By resolving biomolecular structures and probing the dynamics of biomolecular processes, both in vivo and in vitro, researchers have developed a global understanding of how a cell functions.

A limitation of these experimental techniques is the spatio-temporal resolution that they can probe, particularly within the complex and crowded environment of the cell. Molecular dynamics (MD) simulations provide a suitable alternative approach, covering the relevant length and timescales at molecular resolution, albeit over short periods of a cell cycle. Over the past decades, MD has matured into a powerful tool that functions as a computational microscope (Lee et al., 2009; Dror et al., 2012). With the advances in available computer power, including the transition from using central processing units (CPUs) to graphical processing units (GPUs), the complexity and the spatio-temporal scales of MD simulations have increased remarkably. State-of-the-art simulations, containing hundreds of millions of atoms, include dynamic models of a photosynthetic chromatophore vesicle (Singharoy et al., 2019), the nuclear core complex (Mosalaganti et al., 2022), the membranes of a mitochondrion (Pezeshkian et al., 2020), the bacterial cytoplasm (Yu et al., 2016), and a virus particle embedded in a nanoscale aerosol droplet (Dommer et al., 2022). The natural next step is, arguably, the scale of entire cells (Bhat and Balaji, 2020; Khalid and Rouse, 2020; Luthey-Schulten et al., 2022; Thornburg et al., 2022).

Creating a whole-cell model has long been a major goal in computational modeling. A computational cell will help us to resolve a more integral picture of how biomolecular interactions drive cell function since biomolecular processes operate on a hierarchy of interconnected scales. Thus, resolving the full cell function requires a holistic approach. The current state-of-the-art uses static representations of heterogeneous cell-scale structures such as cellPACK (Johnson et al., 2015; Maritan et al., 2022), genome-scale well-stirred reaction models for metabolism (Karr et al., 2012; Macklin et al., 2014; Karr et al., 2015; Breuer et al., 2019; Macklin et al., 2020), or mesoscale kinetic models that attempt to include both structural and chemical states of the cell such as Lattice Microbes (Roberts et al., 2013). These computational techniques, despite granting impressive insights into the complexity of cellular processes, are limited by the spatial resolution they can achieve.

Constructing whole-cell models requires the integration of a large amount of experimental data, i.e., relies on an integrative modeling approach (Bonvin, 2021; Luthey-Schulten, 2021; Gupta et al., 2022). Obtaining such data with high spatial and dynamic detail, particularly in living cells, is very challenging, but exciting progress is being made in elucidating the architecture and stoichiometry of cellular components at unprecedented resolutions (Reading et al., 2017; Ando et al., 2018; Cheng, 2018; Chorev et al., 2018; Christie et al., 2020; Lorent et al., 2020; Narasimhan et al., 2020; Wietrzynski et al., 2020; Štefl et al., 2020), setting the stage for spatially detailed simulations of whole cells. To showcase this possibility, we consider one of the simplest cells known to date: the minimal cell created by the J. Craig Venter Institute (Hutchison et al., 2016), a stripped-down version of a Mycoplasma bacterium. The current strain, named JCVI-syn3A (Syn3A), contains only 493 genes and is still able to replicate independently (Breuer et al., 2019). This cell is particularly amenable to detailed computational modeling because it is of relatively small size (measuring 400 nm in diameter), and its precise composition has largely been resolved (Breuer et al., 2019).

Here we present our ongoing effort to simulate Syn3A using the Martini coarse-grained (CG) force field (Marrink et al., 2022). The Martini force field employs a four-to-one mapping scheme, where up to four heavy atoms and associated hydrogens are represented by one CG bead. This reduction in the number of particles in the system, together with a smoothening of the potential energy surface, speeds up the simulations by about three orders of magnitude, enabling simulations of systems that approach the size of entire cells. In the case of Syn3A, about 550 million CG particles are required, corresponding to more than six billion atoms.

In the remainder of this work, we first describe the set of tools needed to construct a system as complex as an entire cell at the Martini level (the Martini “ecosystem”), including a description of the stepwise procedure to construct a starting model for Syn3A, from building the chromosome and modeling the cytoplasm to the addition of the cell envelope. We end with the prospects of actually simulating this model and discuss the potential future avenues of simulating entire cells. The integrative modeling workflow is schematically depicted in Figure 1.

Modeling cellular environments using a coarse-grained approach requires the use of a force field that incorporates enough detail to represent all biomolecules and their interactions explicitly. In this regard, the Martini force field is an excellent candidate, as demonstrated by the wide range of application studies using Martini over the past 2 decades (Marrink et al., 2022). Additionally, parameters for a large variety of biomolecules are already available, including proteins (de Jong et al., 2013), lipids (Wassenaar et al., 2015), polynucleotides (Uusitalo et al., 2015; Uusitalo et al., 2017), carbohydrates (López et al., 2009; Grünewald and Punt, 2022) and metabolites (Sousa et al., 2021; Alessandri et al., 2022). A curated collection of all parameters is available from the Martini Database (Hilpert et al., 2022).

Accompanying the Martini force field is a collection of tools that compose the software side of the Martini ecosystem (Figure 1). The primary goal of this software suite is to facilitate the construction of topologies and initial coordinates for running CG Martini MD simulations. The Martini ecosystem is currently being extended around a central framework, named Vermouth. Making use of a graph-based description of molecules, Vermouth implements a unified handling of processes frequently used in Martini, such as topology and coordinate generation or resolution transformation, as a lightweight python library (Kroon et al., 2022).

Proteins comprise the bulk of a cell’s organic material, composing approximately 40% of the total intracellular volume (Ellis and Minton, 2003). The number of unique proteins expressed by the cell, i.e., the proteome, can range from a few hundred to several thousand. Consequently, describing realistic cellular environments requires generating topology files for proteins of varying shapes and sizes as well as packing these into a highly crowded solution. Martinize2, which is built on top of Vermouth, is designed for high-throughput generation of Martini topologies and coordinates for proteins from atomistic protein structures. The workflow used in Martinize2 additionally performs quality checks on the atomistic protein structures and alerts to potential problems making it ideal for such a high-complexity use case (Kroon et al., 2022). To generate dense protein solutions on a cellular scale as required for whole-cell modeling, a new software tool, called Bentopy, is currently under development. It uses an efficient collision detection scheme (Howard et al., 2016) to generate random packings of proteins and protein complexes within volumetric constraints. Furthermore, functional annotations of proteins can be integrated into the algorithm, biasing their spatial distribution in the cytosol based on their known biochemical function.

Constructing coordinates and input files for chromosomal DNA presents another challenge for modeling of a whole cell. Even for a comparatively small genome as that of JCVI-syn3A, approaches that rely on reading input coordinates and forward mapping such as used in Martinize2, become too inefficient due to the sheer size of the molecule. Instead the Polyply software, which was originally developed to efficiently setup general polymeric systems, can be used (Grünewald et al., 2022). Within Polyply, a multiresolution graph-based approach is used to efficiently generate polymer topologies, in particular for DNA, only from the sequence. In addition to topologies, system coordinates can be generated using a specialized biased random walk protocol. This tool of the Martini ecosystem has successfully been applied to model dense polymer melts and simple ssDNA viral chromosomes. At the moment, the package is being extended to handle double-stranded nucleic acids, and generate more complex DNA structures such as bacterial chromosomes.

Lastly, modeling lipid membranes has historically been a leading application of the Martini force field (Marrink et al., 2019). Simulating arbitrarily complex membranes of various sizes, geometries, and lateral heterogeneities is facilitated by TS2CG (Pezeshkian et al., 2020). This tool implements a backmapping algorithm that converts triangulated surfaces into CG membrane models. As a result of the method’s high level of control, the curvature-dependent lipid concentrations in both membrane leaflets can be precisely determined by the user. In addition, proteins can be inserted into the membrane together with their characteristic lipid shells [i.e., lipid fingerprints (Corradi et al., 2018)], setting the stage for building cell envelopes.

In the following subsections, we describe the application of the aforementioned tools to construct a proof of principle whole-cell simulation of Syn3A, illustrating how the current Martini ecosystem enables users to study multi-component systems at the mesoscale.

In the wake of a continuous rise in computing power, MD simulations have transitioned from studying idealized representations of biomolecular systems to modeling their full complexity. The culmination of this development would be simulations at the level of entire cells. As a proof of principle that we are ready to meet this challenge, we presented a model of the complete minimal cell JCVI-syn3A, constructed using the Martini ecosystem. The final simulation box comprises more than 560 million CG beads, representing over six billion atoms in the cell (Figure 2).

Before looking at the broader prospects of this endeavor, it is important to discuss a number of limitations of our approach. The current model uses the Martini 2 version of the force field since Martini models for nucleic acids, and other essential cellular components are still under development for the latest Martini 3 release. However, the methods described in this paper can be straightforwardly transferred to the latest version of Martini when validated models become available. With over 800 different bead types and a recalibrated interaction matrix, Martini 3 offers an improved framework for CG MD simulations (Souza et al., 2021). Nevertheless, inherent limitations of Martini, such as an inability to sample protein secondary conformational changes, remain. We do not anticipate that such changes are of primary importance in determining the cellular organization, but details of protein-protein and protein-lipid interactions might be affected. This problem could perhaps be resolved by using Go potentials (Poma et al., 2017; Souza et al., 2019), which are already integrated into Martinize2 and Bentopy.

Even though our in silico cell contains more than 500 unique CG molecules and thereby presumably qualifies as the most complex system simulated to date, it simplifies the composition of various cellular components of Syn3A. Firstly, limited by the availability of Martini models for the metabolites, only a small subset is currently included in the cytosol. Future iterations of our whole-cell model will include Martini models for the complete metabolome, which comprises about 188 different compounds, and are expected to benefit from the ongoing development of dedicated automatic topology builders (Bereau and Kremer, 2015; Potter et al., 2021). Secondly, since AlphaFold2 was used to predict the protein structures of the whole proteome, only monomeric structures were initially available. Essential multimeric proteins like the ribosomes and membrane-embedded transport complexes are either left out or represented by homologous proteins with available experimental crystal structures. In the future, improved protein structure prediction algorithms will be used that also facilitate the modeling of multimeric protein structures. In addition, ongoing progress in the experimental characterization of the Syn3A proteome and lipidome, as well as the characterization of the spatial distributions of membrane proteins, will help further increase our model’s realism. A “living” list of the complete composition of our in silico cell can be found in our GitHub repository (marrink-lab, 2022).

Another issue is the fine-tuning of the amount of interior solvent (both water and ions), together with the lipid balance between the inner and outer leaflet. Previous works on large-scale membrane-enveloped systems (Pezeshkian et al., 2021; Vermaas et al., 2022) have shown that finding this balance is a non-trivial task. Unbalanced systems might experience strong osmotic pressures and membrane (curvature) stress, causing unwanted shape deformations all the way to membrane rupture. As a complicating factor, these effects may only appear after prolonged simulation times. Clearly, dedicated computational resources are required for the simulation of whole cells or cell organelles. The forthcoming generation of supercomputers and simulation software is becoming increasingly efficient, and billion-particle simulations have already been achieved (Jung et al., 2019; Castagna et al., 2020).

An important challenge is reaching timescales long enough to allow meaningful analysis of such large systems. Assuming dedicated computer time on current infrastructure, we anticipate that we can reach timescales of the order of 10–100 µs in the foreseeable future. Although this is typically considered a long enough simulation time for standard system sizes (e.g., a single membrane protein), it is clear that on the scale of an entire cell, we will not be able to equilibrate our system; the generated ensemble of configurations will remain dependent on our starting state. Equilibration will only happen locally, and multiple replicas will need to be generated to obtain statistically relevant data. Note that the 10–100 µs range offers a nice overlap with state-of-the-art experimental techniques. For example, advanced MINFLUX microscopy from the Hell lab enables the tracking of particles as small as 1–2 nm for 100 s of microseconds (Eilers et al., 2018; Schmidt et al., 2021). Besides, Lattice Microbes simulations of the Luthey-Schulten group (Roberts et al., 2013) use time steps of the order of microseconds, which allows for a potential feedback loop between these computational approaches.

Another major limitation is the fact that real cells operate out-of-equilibrium, driven by the import and export of nutrients and an intricate metabolic network of chemical reactions. In our approach, which is based on classical MD, we do not take this into account. We are therefore limited to studying non-reactive processes, i.e., those arising from the physical interactions among the constituents. The current composition of our cell is based on average concentrations of proteins and metabolites and thus reflects a steady-state. Coupling our classic approach with approaches taking into account reactivity, such as the aforementioned Lattice Microbes simulations or other metabolic network models (see below), in principle, could capture the non-equilibrium aspect of real cells.

Keeping these limitations in mind, simulations of the minimal cell with a molecular resolution will make it possible to study a wide range of new aspects. Modelling cellular processes and chemical transformations involves a hierarchy of interconnected scales that cannot be separated without causing artefacts. Behaviour emerging from the interaction of millions of different compounds is easily missed when systems are simplified. One might question to what extent one part of the cell affects another, given the limited timescales likely to be reached. If the various cellular subsystems act independently, one might better simulate those in isolation. To find out, one needs to simulate the complete system in addition to the smaller-scale subsystems. Our whole-cell simulation is only a first step, which will benefit from imminent improvements in high-performance computing to extend these simulations to longer timescales, up to the point where all parts of the cell may influence each other. Currently, the internal organization of the cytosol of Syn3A is primarily a black box. Our model will allow us to observe how proteins inside the cytosol interact with macromolecular structures such as ribosomes and chromosomes. Viewing the cytosol from this perspective, we can observe emerging heterogeneities and viscosity gradients, following in the footsteps of other realistic models of the cytoplasm of various cell types (McGuffee and Elcock, 2010; Yu et al., 2016; Oliveira Bortot et al., 2020). We can expect arising interaction patterns between proteins and metabolites, and probe the possible appearance of biomolecular condensates (Guilhas et al., 2020; Rhine et al., 2020).

A simulation at the level of the entire cell allows us to characterize the extent to which the cell membrane affects (and is affected by) the cellular interior. If we consider a membrane zone with a thickness of 30 nm (∼20 nm of the membrane together with its embedded proteins, plus another 10 nm layer underneath), 40% of the total cell volume is part of this membrane zone. Our simulations will provide detailed insights into the nature and extent of depletion or crowding layers, and into the level of heterogeneity inside this membrane zone, providing information on the extent to which compounds are either enriched or depleted near the cell surface (Nawrocki et al., 2019). A full-cell membrane model might explain why the minimal cell grows on a diet of both saturated and unsaturated fatty acids, but not on a diet of just saturated ones as observed in lipidomics experiments from the Saenz lab (private communication). A related question is why the cell membrane contains such a high percentage of cholesterol (20%–60% dependent on growth medium); this is uncommon for bacterial membranes although generally Mycoplasma do contain some cholesterol for membrane stability.

Of special interest is the potential existence of dynamic highways, i.e., regions in the cell with greater mobility of the constituents, which may arise from crowding effects or liquid-liquid phase separation phenomena, or may be induced by proximity of the cell membrane. Such dynamic highways could be important in regulating transport in an otherwise glassy state of the cytoplasm. For regions of the cell showing particularly interesting behaviour, smaller systems can be extracted with the advanced TS2CG tool and simulated for extended timescales to increase the statistical relevance. Besides passively studying the cellular environment, holistic cell modeling poses the ideal computational sandbox in which we can introduce new components to the cellular environment. For instance, elucidating the non-specific interactions between the cytosol and drug candidates and showing how drug-receptor interactions affect the entire cell instead of just the receptor site.

Using a multiscale modeling approach, we could potentially explore cell dynamics at various stages in its life cycle. Compared to MD simulations, other low-resolution modeling approaches can more broadly explore timescales of several orders of magnitude longer. Integrating other computational models will make it possible to sprout MD simulations in interesting regimes observed with the lower-resolution models. The primary computational method we will focus on integrating into our framework is the whole-cell fully dynamical kinetic model developed by the Luthey-Schulten lab, which accounts for the metabolic pathways governing the cellular processes (Thornburg et al., 2022). By transferring structural information from the kinetic model into our high-resolution model, it will be possible to paint a more detailed picture of the cell’s internal organization and dynamics at specific points of the cell’s life cycle, including during cell fission.

Since most of the tools in the Martini ecosystem are force field agnostic, the workflow can also be applied to generate all-atom whole-cell models. Given the substantial increase in associated computational costs, it might be a wiser approach to only sample smaller subsystems at the all-atom level. These could be straightforwardly obtained from backmapping representative regions taken from the whole-cell CG model. A number of such backmapping tools, optimized for Martini, already exist (Louison et al., 2021; Vickery and Stansfeld, 2021; López et al., 2022).

A final challenge lies in the analysis and interpretation of the complex high-dimensional massive data that will be generated. Clearly, it will be impossible to perform a comprehensive analysis on a whole-cell trajectory, and one needs to focus on specific research questions. However, the trajectories can nowadays be easily shared with the broader community via dedicated open-access repositories such as Zenodo (https://zenodo.org/), allowing others to perform whatever additional analysis they fancy. One can also envision the usage of data reduction schemes to efficiently analyse the whole-cell simulation. One possibility is storing only centers-of-mass movement of the non-aqueous components, which would facilitate the analysis of diffusional behavior, for instance. Another approach would be using a voxel-based method (Bruininks et al., 2021) to dynamically segment the whole-cell model into similarity regions, e.g., membrane periphery or chromosomal region. The system segmentation would allow for efficient quantitative comparison of the cytosolic properties within and between distinct regions of the cell. Furthermore, machine-learning can be invoked to extract interaction patterns and other emergent behavior that might be missed by standard analysis tools (Noé et al., 2020; Wang et al., 2020; Kaptan and Vattulainen, 2022). We foresee that our data sets will generate novel ways of dealing with this unprecedented level of complexity.

In conclusion, we presented a roadmap toward whole-cell MD simulations, illustrated with the construction of the first MD model of an entire cell using our Martini ecosystem. The model represents a next level realized with the computational microscope, providing a complete picture of the cell and making it possible to relate molecular structures and interactions to cellular function directly. In the long term, our computational framework will enable us to study a wide variety of mesoscopic systems, possibly informing the design of fully synthetic cells (Olivi et al., 2021) and modeling cells with more complex internal structures.