1. Introduction

A particularly fascinating problem in science is the description of phenomena involving multiple spatial or temporal scales which together determine the behavior of a system. Many problems in materials science are connected with this structure-property paradigm. The atomic interactions at the microscopic level on the scale of nanometers and femtoseconds fully determine material behavior at the macroscopic scale (on the order of centimeters and milliseconds and beyond). The idea of performing material simulations across several characteristic length and time scales has therefore obvious appeal not only for fundamental research but also as a tool for technological innovation, see e.g., the discussions in ^[1,2,3].

For example, in engineering, understanding the key mechanisms of failure in materials on various length-and time scales is a prerequisite for the design of new materials with desired superior properties such as high strength or toughness. As "material" we understand all condensed forms of matter, i.e., a state of matter, where the electrons are not separated from the atoms (a state which is called "plasma") and thus either form hard matter, i.e., crystalline solids, or soft matter, i.e., biological macromolecules and cells. Most of biological cells are about $1$–$100$ microns in size, and all cells use a plasma membrane to separate and protect their chemical components from the outside environment. The plasma membrane acts as a barrier that separates the interior of the cell from its exterior, preventing the contents of the cell from mixing with the surrounding medium. It also allows certain molecules to pass inward—across the membrane through channel—and waste products to pass out. A plasma membrane is very flexible and deforms without tearing when a cell grows and changes shape.

The interior of the cell is made up of a liquid phase (cytosol), a nucleus, the cytoskeleton consisting of many networks of microtubules, actin and intermediate filaments, organelles of different sizes and shapes, and other proteins. The resistance of single cells to elastic deformation, as quantified by a huge range of the effective elastic modulus as displayed in Figure 1, is orders of magnitude smaller than that of typical engineering materials such as metals and ceramics. This is the reason why cells and polymers are "soft" and governed in their behavior mostly by the interplay of entropy and energy. If a cell's plasma membrane is penetrated, it neither simply collapses, nor does it remain torn; instead, it has the capability of resealing as has been shown in several single cell experiments ^[4]. This transient permeability of a cell's plasma membrane can also be achieved by exposing cells to shock waves, as I have demonstrated in a recent publication where a new experimental set-up for exposing U87 glioblastoma tumor cells to laser-induced shock waves was introdced ^[5]. A shock wave treatment of cells might make possible new interesting biological and medical applications, e.g., in gene therapy, where the permeabilization of the membrane aims at delivering therapeutic molecules, genes or other genetic material into the cells' cytoplasm ^{[6,7,8,9,10,11,12,13,14,15,16]}.

This review article is organized as follows: In Section 2, I provide a review of the basics of shock wave physics. In Section 3, I discuss different aspects of multiscale modeling of materials. Section 4 introduces the concept of coarse-graining which is often (but not solely) used for modeling soft matter systems on the mesoscopic scale. Then, in Section 5, I discuss several recent applications that illustrate how different modeling approaches may be usefully combined in detailed computational studies. Finally, in Section 6, a few general conclusions are drawn and several considerations are offered.

---

2. What is Shock Wave Physics?

The theoretical foundations of shock wave physics originated from 17th-century studies of classical acoustics and 18th-century aeroballistics. With the emergence of high-speed photography in late 19th century and further progress in experimental visualization techniques in 20th century, the study of shock waves in different states of matter has gradually emerged from a very small and unnoticed branch of physics to a major complex and interdisciplinary science ^[18,19].

Shock waves are in essence discontinuous, entropy-increasing, rapid mechanical phenomena that have been observed and studied in the laboratory and in nature, in microscopic as well as in macroscopic dimensions and in all states of matter: gaseous ^[20,21], liquid ^[22,23,24], solid ^[25,26,27], plasma ^[28,29], and even in Bose-Einstein condensates ^[30,31,32]. Terrestrial examples ^[33] of naturally occurring discontinuities are high-energy events such as volcanic explosions, thunder, meteorite impacts, sea-and earthquakes or tsunamis, while in outer space ^[34,35] they encompass plasma shock waves induced by solar wind, supernovae explosions, implosions of white dwarfs, comet and asteroid impacts, and stellar or galactic jets. This variety of phenomena and in particular the possible applications of shock waves in different areas such as biology, chemistry, physics, medicine and even engineering renders the scientific study of shock waves a rather fascinating subject.

---

2.1. Definition of Shock Waves

A shock wave is a mechanical wave characterized by an area or sheet of discontinuity in which, within a narrow region, thermodynamic quantities such as pressure $p$, density $\rho$, particle velocity $\vec{v}$ or temperature $T$ change abruptly. While a theory of shock waves with an infinitesimal jump condition can be developed mathematically in the framework of approximative continuum models of condensed matter, i.e., based on the hydrodynamic equations of ideal fluids or gases, in physical applications, shock wave theory is somewhat limited, because the width of a genuine shock wave is always of finite size and because of the actually discrete atomic nature of condensed matter. It turns out that shock waves are in essence small regions were non-adiabatic, i.e., irreversible energy dissipation occurs and for which the shock wave is a mathematical idealization. The width of the shock wave, i.e., the size of the dissipative region, establishes itself according to the conservation laws of continuum theory: the conservation of mass, momentum and energy. In very strong, high-amplitude shock waves, their width becomes so small that they are practically indistinguishable from the mathematical idealization of an infinitesimally small perturbation with jump conditions for thermodynamic variables, see Figure 2.

In solids, physical shock waves are mechanical waves of finite amplitudes and arise when condensed matter is subjected to a rapid compression. Shock waves can be defined by several major distinctive properties as follows:

• the formation of a steep wave front with a sudden, irreversible change of thermodynamic quantities;

• a pressure-dependent, supersonic velocity of propagation;

• the creation of entropy;

• non-linear superposition properties (for reflection and interaction).

Sometimes another criterion for a shock wave is listed in shock wave literature which is an extremely short rise time of the pressure within tens of nanoseconds. However, this is not a defining criterion of shock waves, because no generally agreed upon definition exists of what "extremely short" really is supposed to mean—it is a term of rather subjective nature. A nanosecond rise time certainly applies for lithotripter shock waves which are used in medicine as a non-invasive technique to comminute calculi ^[36] and for other therapeutic purposes ^[37]; however, in the geologic realm for example, meteorite impacts, explosive volcanic eruptions and earthquakes can provoke drastic and irreversible changes within seconds which is nine orders of magnitude slower.

One thing that distinguishes a shock wave from an ordinary sound wave is that the initial disturbance in the medium that causes a shock wave is always traveling at a velocity greater than the phase velocity of sound in the medium. The electromagnetic analog of this is called Cherenkov radiation ^[38], which is caused by a charged particle traveling through a dielectric medium at a velocity faster than the speed of light in that medium. Because a shock wave moves faster than the speed of sound, the medium ahead of the shock front cannot respond until the shock strikes, and so the shock wave falling upon the particles of matter initially at rest, is a supersonic phenomenon. In a steady flow with constant velocity $\vec{v}$ the speed at which a perturbation travels relative to a reference frame in which the laboratory is at rest, is composed of two parts, see Figure 3. On the one hand the perturbation travels with velocity $\vec{v}_s=\vec{c}$ relative to the fluid into a direction $\vec{n}$. On the other hand, the source of disturbance moves at the same time with the velocity $\vec{v}$ of the fluid flow. Figure 4 exhibits one of the first ever published schlieren picture of the bow shock wave generated by a projectile moving at supersonic speed. The two vertical lines visible in Figure 4 are wires with a spark gap inserted for triggering the spark light source. The projectile had a speed of $v_0\approx 530\, m/s$ which resulted in a Mach cone angle of $\alpha\approx 40^{\mathrm{\circ}}$ (cf. Figure 5b).

For simplicity, we assume here constant flow velocity $\vec{v}$ and consider a perturbation in the fluid at a fixed point $O$. This perturbation propagates with velocity $\vec{v} + c\vec{n}$ from $O$. We can see in Figure 5 that the velocity depends on the direction of $\vec{n}$. Vector $\vec{n}$ is a radial unit vector starting from point $O$ with respect to the moving fluid flow. The possible values of the velocity lie on the surface of a sphere according to Figure 5. If $v_s < c$, the vectors $\vec{v} + c\vec{n}$ can attain any direction in space, cf. Figure 5a. Hence, a flow that propagates with subsonic speed will affect the whole domain of the fluid. For $v_s > c$, the velocities $\vec{v} + c\vec{n}$ lie within a cone the tip of which is located in point $O$, cf. Figure 5b. This cone is the envelope of the sphere, so we find for the total cone angle, the so-called Mach angle $\alpha$ the following relation:

$ \sin\alpha = \frac{c}{v}. $

(1)

The ratio

$ M_a = \frac{v}{c} $

(2)

is called Mach number. For $M_a > 1$, the flow is supersonic and for $M_a < 1$ it is subsonic, respectively.

Hence, we find that a perturbation travels only with supersonic speed in the direction of flow within the boundaries of a cone with angle $2\alpha$. Angle $\alpha$ gets smaller with increasing flow speed $v$, so the perturbation in a flow has no effect in the region outside that cone. The boundary layer of the region which can be influenced by a perturbation is called Mach's cone, characteristic lines or simply characteristics. The parts c) and d) in Figure 5 exhibit the situation in a fluid at rest through which a perturbation moves with velocity $\vec{v}$.

Another difference between a shock wave and an ordinary sound wave is that the entropy is increased in a shock wav—hence, in contrast to ordinary sound waves, shock waves constitute an irreversible process. In engineering applications, shock waves produced in air by an explosion and radiating outward from its center are termed blast waves, because they cause a strong wind, while the term shock wave is preferred for such waves occurring in water or the ground, because here the effect is like that of a sudden impact. Impact experiments are often performed to study the fracture and failure behavior of solids or in an attempt to determine a shock equation of state, describing a material under extreme conditions of high pressure and/or temperature. Along with phenomenological constitutive equations, this allows for numerically analyzing a material in a shocked state, e.g., in a finite element analysis, cf. Figure 6.

Upon hydrostatic compression, a solid shrinks in volume and the pressure increases ^[41]. The hydrostatic curve is smooth as a function of stress and volume. If the solid is compressed uniaxially, stress and pressure also increase until the maximal shear stress of the material is reached. Then it will yield and change from the uniaxially compressed state to hydrostatic compression. Many solid materials under moderate shock loading exhibit this two-wave structure where the first one is essentially an elastic wave followed by a slower plastic deformation wave front, which can be seen in Figure 6, too. In the $p-V$-diagram this transition from elastic to plastic behavior is visible through a cusp, called the Hugoniot elastic limit (HEL), and there are a considerable amount of experimental data which illustrate these elastic–plastic–shock phenomena (EPSP) with an elastic precursor wave ^[42].

Very often, impact experiments are used as test cases for the predictions obtained from numerical multiscale models—such as the ones presented in Figure 6—which simulate a material specimen that can be macroscopic in size and take several micro-structural features of the specific solid into account. Another important use of impact experiments in the laboratory is to study the crack initiation and propagation until ultimate failure occurs in a set-up that allows for very reproduciblex shock wave generation. One of these set-ups is the so-called edge-on-impact test, see Figure 7, where an impactor strikes the edge of a material specimen at high speed, initiating a well-defined plane shock wave through the material that ultimately leads to fracture and destruction.

Another important application of research in shock wave physics is the investigation of hypervelocity impact phenomena (with relative speeds of objects beyond $\sim2$ km/s) which occur in impact events in space. Here, scientists try to understand via novel computer simulations and impact experiments on a laboratory scale the physics behind the formation of space debris clouds in the low Earth orbit (LEO) ^[44]. Such clouds form when small debris objects from former space craft and satellite missions impact a satellite in its orbit at hypervelocity speed. These types of investigations also contribute to the assessment of the socio-economic consequences that are involved in the risks in space from orbiting debris.

Finally, also phase transitions in solids can be induced using shock waves and have been studied experimentally and numerically using classical molecular dynamics simulations in a variety of materials ^{[26,45,46,47,48]}. The world record for the largest ever MD simulation of shock waves in a solid is currently still held by Germann and Kadau ^[49] from Lawrence Livermore National Laboratories, who simulated a simple cubic crystal consisting of $1$ trillion ($10^{12}$) atoms, treated as classical Lennard-Jones particles. This corresponds to a chunk of matter of edge length of about $2.5\, \mu m$.

---

3. What is Multiscale Modeling?

In recent years, we have seen a very large increase in computing power and—at the same time— progress in developing efficient algorithms. As a result, computer simulations have become very popular and are used today in almost all areas of research ^{[50,51,52,53,54]}. Fast parallelized algorithms allow for the solution of non-linear many body problems directly, without any mathematical approximations involved. Usually, in analytic theory, practically all non-trivial problems are way too complex to be solved with just pencil and paper. In this sense, computer simulations can provide a link between experiment and analytic theory, thus allowing to test theories. But simulations can also be used as an exploratory research tool under physical conditions not feasible in real experiments in a laboratory. Thus, doing science with the help of computer systems has established a new, interdisciplinary research approach which is often referred to as "Computational Physics" or "Computational Materials Science". This approach combines knowledge and methods of different research areas such as chemistry, biology, physics, mathematics, engineering and even medicine and allows for performing multiscale and multi-disciplinary simulations in situations relevant for real-world applications.

Nowadays, with fast computational tools at hand, simulations have become feasible that are of practical interest in engineering sciences for product design and testing. Many materials used in industrial manufacturing are very heterogeneous and can be characterized by different kind of defects, interfaces, and other microstructural features. As an example, in Figure 8 we display a Scanning Electron Microscope (SEM) photograph of the fracture surface of Aluminum Oxide (Al$_2$O$_3$) after planar impact load. In general, inorganic crystalline materials have structural features such as grain boundaries between crystals which are millimeters to micrometers in size, while on the atomic scale dislocations and defects such as vacancies occur. Hence, these structures have to be studied from a hierarchical perspective.

For example, simulations in material physics and chemistry often focus on the exploration of lattice and defect dynamics on the atomic scale using MD and Monte Carlo (MC) methods based on generalized force-fields (physical potentials), deduced from solving the non-relativistic Schrödinger equation for a limited number of atoms ^[55,56,57]. Figure 9 displays several methods commonly employed for atomic scale simulations along with their estimated maximum system size and the typical time scales that can be treated. The highest precision and transferability of methods is achieved with self-consistent first principles—so-called ab-initio—calculations which try to avoid completely any empirical fitting parameters.

Self-Consistent Field Theory (SCFT) is in essence based on the Hartree-Fock (HF) method which itself is a mean-field approximation (MFA), using MFA for energy calculations always leads to values which are larger than the exact energy values. Another approximation in HF calculations stems from the need to find an analytic expression for the wave function, such functionals are only known exactly for very few one-electron systems and therefore, usually some approximate functions are used instead. Normally, as a simplification, one uses as basis wave functions either plain waves, Slater type orbitals (STO) $\sim \exp(-ax)$, or Gaussian type orbitals (GTO) $\sim \exp(-ax^2)$. The term "orbital" is just another name for denoting a one-electron wave function. Correlations, i.e., interactions between the electrons are treated within Møller-Plesett perturbation theory (MPn), where n is the order of correction, Configuration Interaction (CI), Coupled Cluster theory (CP), or other methods. As a whole, these methods are referred to as correlated calculations.

An alternative method to SCE for calculating ground state energies of molecules is Density Functional Theory (DFT). In DFT, the total energy of a system is not derived from a wave function but in terms of an approximate Hamiltonian and thus from a Local Density Approximation (LDA). DFT uses GTO potentials or plane waves as basis sets.

Tight Binding (TB) is a so-called semi-empirical method (see e.g., Porezag et al. ^[58] or Pettifort ^[59]) which uses an approximation of the Hamiltonian ${\cal H}$ by approximating or neglecting several terms (called Slater-Koster approximation), but re-parameterizing other parts of ${\cal H}$ in a way so as to yield the best possible agreement with experiments or ab initio simulations. In the simplest TB version, the Coulomb repulsion between electrons is neglected. Thus, in this approximation, there exists no correlation problem, but there is also no self-consistent procedure ^[60].

In classical MD methods, one integrates the classical Newtonian equations of motion for a large number of particles, respectively atoms, which are fully treated as classical point particles or as classical particles with tare volume. That is, MD in essence solved the classical $N$-body problem. A separate consideration of the motion of the electrons, as in the ab-initio methods, is not done. Interactions between particles are usually described via empirical, sometimes generic potentials which may be fitted to the results of experiments or of quantum mechanical calculations. Due to the increasing computing power of modern hardware, many-particle MD simulations taking into account the degrees of freedom of several billion atoms are nowadays feasible ^[54,61]. Molecular dynamics simulations of this kind using generic potential models for a solid enhanced our understanding of the basic processes that govern failure and crack behavior, such as the dynamical instability of crack tips ^[62,63], the limiting speed of crack propagation ^[64], the dynamics of dislocations ^[65], and the universal features of energy dissipation in fracture ^[66].

In contrast to the natural sciences, simulations of materials in the area of mechanical engineering typically focus on large-scale problems, where the $N$-body problem is often homogenized with approximate continuum methods, such as Smoothed Particle Hydrodynamics (SPH) or the Finite Element Method (FEM), utilizing additional, purely empirical constitutive relations ^[67]. The SPH method is simply a re-formulation of the hydrodynamics equations for (pseudo-)particles ^[68]. This methods suffers from certain important deficiencies such as convergence problems. Finally, the FEM is simply a continuum approximation of material structures, where the structural details are represented by a topologically connected mesh, where the mesh nodes are the integration points of the continuum conservation equations. Usually, this method reaches its size limitations very quickly and only meshes with a couple of million integration nodes for macroscopic structures can be simulated on large computer systems. One disadvantage of this method is the use of a mesh which is very bad for simulating large deformations in structures and which may lead to prohibitively small time steps in the numerical procedure. Also, modeling failure and fracture with this method is very awkward and leads to the introduction of artificial, unphysical corrections such as artificial viscosity in an attempt to stabilize the numerical procedure.

Generally speaking, in materials with small grain size the diffusion distances are small. Hence, processes governed by diffusion, such as sintering, are facilitated and can occur at lower temperatures than would otherwise be possible. Being able to predict the properties and performance under load of such materials is a major issue in current computational material research and for industrial product design. Unfortunately, the complexity of structural hierarchies in solids on different scales does not allow for selecting one single computational model or physical theory which allows for explaining all observed phenomena. Therefore, the different kind of microstructural features of different important classes of materials such as metals, ceramics, or materials pertaining to soft matter (glasses or polymers), have to be considered in different models, cf. Figure 10.

---

3.1. Hierarchical Length and Time Scales

Solids comprise a span of length scales of roughly $10$ to $12$ orders of magnitude and the methods of classical Newtonian physics are sufficient to describe many of the occurring phenomena, cf. Figure 11. However, classical MC or MD methods lose their validity at length scales comparable to the typical size of atoms ($\approx 10^{-10}$ m). In classical theories atoms are treated as point particles or spheres. In principle, the non-relativistic, time-dependent Schrödinger equation describes the properties of molecular systems with high accuracy, but anything more complex than the equilibrium state of a few atoms cannot be handled at this level. Therefore, approximations have to be introduced. The approximation are the more severe, the more complex a system is and the longer the time scale of the investigated processes are. At a certain point, the ab-initio approach has to be abandoned completely and replaced by empirical parameterizations of the used model. Hence, depending on the kind of problem that one investigates and depending on the desired accuracy with which specific structural features of the considered are to be resolved, one has the choice between many different models for different length and time scales. Unfortunately, there is no simple "hierarchy" that is connected with a natural length scale according to which the great diversity of simulation methods could be sorted out. For example, Monte Carlo lattice methods can be applied at the femtoscale of Quantumchromodynamics ($10^{-15}$ m) ^[69], at the Å ngstrømscale ($10^{-10}$ m) of solid state crystal lattices ^[70], or at the micrometerscale ($10^{-6}$ m), simulating grain growth processes of polycrystal solid states ^[71].

The salient hierarchical structural features of materials have to be taken into account when developing mathematical and numerical models. Here, usually one of two possible strategies is pursued: In a "sequential modeling approach" one attempts to compile together a hierarchy of computational approaches in which large-scale models use the coarse-grained representations with information obtained from more detailed, smaller-scale models ("bottom-up" vs. "top-down" approach), see Figure 11 and also compare Figure 12 in Section 4.

This sequential modeling technique has been used effectively for systems in which the different scales are weakly coupled. The large majority of multiscale simulations are based on a sequential approach. Examples of this approach are abundant in literature and it is not my intention to comprehensively review the many publications in this field. The sequential approach includes practically all MD simulations the potentials of which are derived from ab-initio calculations, including classical coarse-grained simulations of macromolecules ^[72,73,74]. Due to the fractal nature of polymers, they exhibit self-similar structures on various length scales. This leads to unique and universal scaling properties which can be used to check the validity and quality of research codes. This "Scale-Hopping in Computer Simulations of Polymers" ^[75] is typical for coarse-grained simulations of soft matter systems.

The second strategy used in multiscale modeling and simulation is the "concurrent" or "parallel approach". Here, one links appropriate methods on each scale together in a combined model, where the different scales of the system are considered concurrently and often communicate with some type of hand-shaking procedure ^[76]. This approach is necessary for systems where the properties at each scale strongly depend on what happens at the other scales, for example dislocations, grain boundary structure, or dynamic crack propagation in polycrystalline materials ^[77,78].

---

4. What is Coarse-Graining?

The idea of using coarse-grained (CG) models for the description of molecular structures was originally introduced by Warshel and Levitt ^[79] in a 1976 pioneering paper which dealt with the structure of globular proteins. Since then, CG models have found their way into polymer physics as so-called bead-spring models, taking advantage of universal scaling laws of long polymer chains which are due to their fractal nature ^{[80,81,82,83]}, as well as into chemistry, geophysics, engineering and other areas of computational research ^[43,84].

CG models of macromolecules provide a route to explore biomolecular systems on larger length and time scales while still resolving important physical aspects of the molecular topology, for example, the specific lipid bilayer structure of biological membranes ^{[3,72,83,85,86,87,88,89]}. They constitute a class of mesoscale models, in which many atoms or groups of atoms are treated by grouping them together into new particles which act as individual interaction sites usually connected by entropic springs, see Figures 12 and 13.

Such models have become popular ^[92,93] because they remain particle-based while at the same time reducing the computational cost. Another advantage of CG methods is that the groups of particles can be constructed with various levels of detail, thus permitting, e.g., the study of biological membranes on multiple length scales. A "bottom-up" CG model is a model of a particular system that is constructed on the basis of a more detailed model for the same system as indicated in Figure 12. In principle, the high-resolution, all-atom model may be based on atomistic data deduced from atomistic structure calculations. In contrast, "top-down" models do not rely upon or directly relate to a more detailed model for a particular system. Instead, they are usually related to the full complexity of the real experimental system by addressing observables on length scales that are accessible to the CG model. Often, these observables are thermodynamic averaged quantities such as pressure, temperature, stress and strains or forces accessible by direct experimental measurement. Figure 12 illustrates schematically these two major approaches to coarse-graining.

CG simulations involve less computational effort than atomistic simulations, because the number of interacting particles is reduced. Consequently, CG simulations are able to investigate much larger length-and time-scales than is possible in all-atom approaches, let alone in quantum chemical calculations ^[72,94,95]. CG models may simply remove certain degrees of freedom (e.g., vibrational modes between two atoms) or they may in fact simplify the two atoms completely via a single particle representation. The ends to which systems may be coarse-grained is simply bound by the accuracy in the dynamics and structural properties one wishes to replicate. This area of research is still basically in its infancy, and although it is commonly used in biological modeling, chemistry and polymer physics, the analytic theory behind it is still poorly developed.

Recent literature dealing with multiscale methods in polymer physics includes the edited volume Multiscale Molecular Methods in Applied Chemistry ^[74] and the treatment of multiscale computational methods in the monograph Computational Multiscale Modeling of Fluids and Solids—Theory and Applications ^[3]. Additionally, very recent publications on this subject can be found on the website of the special issue Computational Multiscale Modeling and Simulation in Materials Science of the Journal Materials ^[96].

---

4.1. Coarse-Graining of Soft Matter: Polymers and Biomacromolecules

The research area of CG simulations of polymers and biological macromolecular structures has seen major progress over the past $30$ years, which is mostly due to the emergence of more physically based modeling approaches in the biological sciences, which allowed for the exploration of soft biomaterials, enhanced understanding of structures and diseases and which lead to the development of new medical treatments and diagnostic methods ^{[97,98,99,100,101]}. The major challenge that material scientists are facing in this research area is the extremely broad spectrum of length and time scales in the dynamics of biological macromolecules. Because of their large length, macromolecules can exist in many possible conformations that contribute to the total free energy of the system. In essence, the structural properties of macromolecules are determined by an interplay between entropic and energetic contributions. The hydrodynamic interactions of macromolecules with their surrounding solvents, covalent interactions, intra-and intermolecular interactions and—particularly in biological macromolecules—the Coulomb interactions between charged monomers along with hydrogen bonds add up with entropic forces to build a very complex system of interacting constituents. Because of their large range of characteristic time and length scales, macromolecules are used in many technological applications. The theoretical description of such systems has to take into account some of their hierarchical salient features to ensure equilibration on all relevant length and time scales and to efficiently sample their complex potential energy hypersurface.

CG models—as shown in Figures 12 and 13—are used because atomistic models of long polymer chains are usually intractable for time scales beyond nanoseconds, maybe microseconds on the largest computer systems; however, these larger time scales are important for many physical phenomena (for example diffusion processes) and for comparison with real experiments. Macromolecules, due to their large molecular weight exhibit fractal properties. Hence, numerical models of this type molecules are very useful for the investigation of fundamental scaling properties. In fractal systems, there exists some property $P$ of a system (with polymers, e.g., the radius of gyration $R_g$) which obeys a relation $P\propto N^k$, where $N\in\mathbb{N}$ is the size of the system, i.e., the number of monomers and $k\in\mathbb{Q}$ is the fractal dimension which is of the form $k=a/b$ with $a\neq b$, $a\in\mathbb{N}$ and $b\in\mathbb{N}$. In order to extract the universal scaling properties of polymers with CG numerical models it is sufficient to ensure the following features:

• connectivity of monomers in a chain;

• The topological constraints, e.g., the impenetrability of chain bonds;

• The flexibility or stiffness of monomer bonds.

CG models of lipids and proteins have allowed for simulating lipoprotein assemblies and protein-lipid complexes for several microseconds using just a couple of coarse-grained beads for each amino acid of the protein. Without coarse-graining, such simulation would not have become possible.

---

5. Research Examples in Polymer Physics and Computational Biology

Polymers are encountered in a large variety of chemical achievements, from biopolymers like Desoxyribonucleic Acid (DNA), cellulose, proteins, and actin filaments, to synthetized polymers like polyethylene (PE), polybutadiene (PB), polystyrene (PS), polymethylmethacrylate (PMMA)—which is used in the production of compact discs—which are all of great industrial interest. They are employed under numerous forms, in solutions, as gels, rubbers, synthetic fabrics, molded pieces, and so on.

---

5.1. The Ideal Polymer Chain Model

Here, I limit myself to the description of main chain flexible polymers, whose prototypical monomer is polyethylene [-CH₂-]n. This monomer is repeated many times in polymerization reactions to form what is called a polymer, i.e., a macromolecule that consists of many repetitions of one single monomer unit. Other standard examples are provided in Figure 14.

Every single one of these polymers can be represented as a series of bonds linking successive monomers; for the simple case of polyethylene (PE), the angular variation between neighboring carbon atoms is given by $\Theta_{n-1}=0$, $\Delta\varphi_n=[0, \pm 120^{\circ}]$, see Figure 15.

The so-called trans configuration of polyethylene is characterized by $\Delta\varphi_n=0$ and $\Delta\varphi_n=\pm 120^{\circ}$ denotes gauche configurations. Gauche and trans configurations differ from each other energetically. The geometric picture shown in Figure 15 for polyethylene can be extended to any type of homopolymer.

Assuming the successive values of $\Delta\varphi_n$ to be random, the polymer chain appears, when observed at some scale much larger than the typical average bond length $a$ of monomers, as a statistical coil. The persistence length $L_p$ is then a measure for the degree of stiffness in such a polymer coil. $L_p$ denotes the length scale on which the chain looks like a sequence of rigid collinear segments in trans configurations. Let $\Delta\epsilon = \epsilon_g -\epsilon_t > 0$ be the difference in energy between the gauche configuration and the trans configuration. Then, the persistence length is given by

$ L_p = a\exp\left(\frac{\Delta\epsilon}{k_BT}\right), $

(3)

where $T$ denotes temperature and $k_B$ is the Boltzmann constant. The concept of $L_p$ is based on the idea that the correlation between the orientations of the local tangent $t(s)$ decays with the distance along the filament contour according to $\langle t(s)\, t(s^{\prime})\rangle = \exp\left(-\lvert s-s^{\prime}\rvert/L_p\right)$. The total chemical length of the chain is given by $L=Na$, where $N$ is the degree of polymerization, i.e., the number of repeat units (or monomers) in the polymer. $L$ has to be large compared with $L_p$ for the chain to be considered as a statistical coil, i.e., $L_p/L\ll 1$. The degree of polymerization $N$, i.e., the number of repeating monomer units in a polymer is sometimes also called the number of segments which is not to be confused with the bond vectors $\vec{a}_i$ with the convention, that index $i$ refers to the position vector connecting segment number $i-1$ with segment number $i$. Thus, the counting of bond vectors starts with $i=1$ and ends with $N-1$. In a polymer model, one often just refers to the number $N$ as the number of monomers, or—in computer simulations—the number of particles.

An example of polymers with vastly different persistence lengths is displayed in Figure 16 which shows a classification of polymer stiffness of the three major polymer components of the cytoskeleton of eukariotic cells, namely microtubules (stiff rods), actin filaments (semiflexible polymers) and intermediate filaments (flexible polymers).

Phenomenologically, the persistence length can be thought of as a characteristic length of the thermal fluctuations of an elastic rigid rod. Let us assume, for the sake of simplicity, that the cross-section of the rod is isotropic. Then, the free energy per unit length of the rod can be written as:

$ F=\frac{1}{2}\kappa r_c^{-2}, $

(4)

where $r_c$ is the radius of curvature. We have $r_c^{-1}=\pm\lvert\partial^2\vec{u}/\partial s^2\rvert =\partial\omega/\partial s$, where $\vec{u}$ is the small displacement of the rod with respect to the unstrained straight rod, $s$ is the length measured along the rod, and $\omega$ is the angle of deflection of the rod. lt is possible to show, by Fourier transforming the free energy $F=\smallint_0^L\kappa r_c^2\, ds$ of a rod of length $L$, then applying the theorem of equipartition of energy, and finally summing over all modes, that

$ \langle\omega^2\rangle=2\frac{k_BT}{\kappa}L. $

(5)

Taking $\langle\cos\omega\rangle\approx 1 -Lk_BT/\kappa$ equal to zero, which is the condition at which the elements of the rod at $s = 0$ and $s = L$ are decorrelated, one gets

$ L_p=\frac{\kappa}{k_B T}. $

(6)

Thus, for $L < L_p$, the density of the free energy density can be written as

$ F=\frac{1}{2}k_BTL_pr_c^{-2}. $

(7)

We also mention the concept of a persistence time $\tau_p$, i.e., a time characteristic of a change of conformation. Let $\Delta E$ be the activation energy necessary for the trans-to-gauche isomerization. Then we have

$ \tau_p =\tau_0\exp\left(\frac{\Delta E}{k_B T}\right). $

(8)

The polymer chain appears as flexible during any observation time $\tau > \tau_p$. The characteristic time $\tau_0$ is related to the thermal vibrations of the chain, and is of the order of $10^{-12}\ $s.

---

5.2. Single-Chain Conformations

For a polymer chain there are three rotational isomeric states per bond: $t$, $g^+$, and $g^-$. Hence, a polymer chain can attain $3N$ different configurations. Among those, some are special, as follows:

• Helix Configurations. In a helix configuration, all monomers are repeated along a helical elongation of the chain. Helical configurations have minimal internal energy and are observed in crystalline polymers. For polyethylene, this is the all-trans state: The chain exhibits a zigzag structure along its elongation, rather than a true 3D helical structure. For polytetrafluoroethylene (PTFE), it is a true helix, not far from the all-trans state, but differing from it by a superimposed twist along the elongation, caused by the repulsive interactions between neighboring fluorines. For polyoxymethylene (POM), there are several minimal helical energy states, obtained by twisting the all-gauche configurations in slightly different ways.

• Coil Configurations. Here, we will occupy us only with molecules that are dynamically flexible, i.e., whose chemical length $L$ is large compared to $L_p$, and whose observation is made on durations much larger than $\tau_p$. These molecules run through the 3N configurational states. In such a case, a number of properties are not dependent on the local scales, and the chemical features can be ignored. On the other hand, all mesoscopic properties depend on $N$ according to scaling laws, of the type

$ R =aN^{\nu}, $

(9)

where $R$ is some quantity that describes the spatial extension of a molecule.

---

5.3. The Ideal (Gaussian) Chain Model

Assuming that the $N$ elementary segments of a polymer chain have no interaction and are independent one from the other. Such a chain is called ideal. The vector that points from one end of the chain to the other one is called end-to-end vector $R_e$ and can be written as

$ \vec{R}_e=\sum\limits_{i=1}^{N}\vec{a}_i. $

(10)

Taking the square of $\vec{R}_e$, and averaging over all possible configurations yields the mean squared end-to-end vector:

$ \langle\vec{R}_e\rangle=\sum\limits_{i=1}^{N}\vec{a}_i{}^2 + 2\sum\limits_{i\neq j}^{N}\langle\vec{a}_i\vec{a}_j\rangle = \sum\limits_{i=1}^{N}\vec{a}_i{}^2 + 2\, \lvert\vec{a}_i\rvert\, \lvert\vec{a}_j\rvert \underbrace{\sum\limits_{i\neq j}^{N}\langle\cos(\measuredangle(\vec{a}_i, \vec{a}_j)\rangle}_{=0} = Na^2, $

(11)

where we have used the (ideal) condition that segments $i$ and $j$ and thus, $\vec{a}_i$ and $\vec{a}_j$ are statistically independent. Hence, the average separation of the ends, which can be viewed as the average size of the chain, is

$ R_e=R_0=\langle\vec{R}_e{}^2\rangle^{1/2}=aN^{1/2}, $

(12)

and exponent $\nu=1/2$. Subscript $0$ in Eq. (12) indicates ideal chain behavior. For non-ideal chains, the expression $R_e$ is used for the average size. Note, that here, $R_e\neq\lvert\vec{R}_e\rvert$ as one might expect from the usual use of the same variable for denoting a vector $\vec{v}$ and its absolute length $\lvert \vec{v}\rvert$.

The scaling law of Eq. (12) remains unchanged when introducing interactions between the different monomers along the same chain (generally called bonded interactions, in contrast to unbonded interactions, e.g., between monomers of different chains). One can show as a consequence of the central limit theorem that the probability $P_N(\vec{r})$ for the $N$-link chain to terminate at $\vec{r}=\vec{R}_e$ is Gaussian distributed and given by:

$ P_N(\vec{r})=\left(\frac{3}{2\pi R_0{}^2}\right)^{3/2}\exp\left(-\frac{3}{2}\frac{r^2}{R_0{}^2}\right). $

(13)

---

5.4. Crossover-Scaling of Linear, Semiflexible Polymers

In this section, we discuss briefly an example of CG modeling of semiflexible polymers. Semiflexibily, i.e., stiffness of a polymer chain is a hallmark of all biological macromolecules, which is why it has to be introduced into modeling schemes for biological macromolecular systems. This can be done with a bending potential that provides an energetic penalty in case the bond angle (defined by the angle between consecutive bond vectors connecting the monomers along a polymer chain) deviates from the preset angle.

In contrast to theoretical modeling, only very few simulation studies using mostly detailed atomistic models have been performed on the behavior of semiflexible polymer systems ^{[102,103,104,105]}. In particular the results of atomistic simulations have shown that the mean square particle displacements and the scattering functions deviate significantly from what is expected under ideal (Gaussian) conditions ^[102].

The Kratky-Porod chain model (or worm-like chain model, WLC) ^[106] is a simple description of semiflexible, inextensible polymers with spatial fluctuations that are governed by their bending energy $U_{\text{bend}}$

$ U_{\text{Bend}} = \frac{\kappa}{2}\int_0^L ds\left(\frac{\partial^2 {\vec{r}}}{\partial s^2}\right)^2 $

(14)

of the inextensible chain of length $L$ depends on the local curvature of the chain contour $s$, where $\vec{r}(s)$ is the position vector of a mass point (a monomer) on the chain and $\kappa$ is a constant ^[107].

A dynamic equation for the WLC model was formulated by Harris and Hearst by applying Hamilton's principle with the constraint that the second moment of the total chain length be fixed. The resulting expressions for the bending ${\vec F}_{\text{Bend}}$ and tension forces ${\vec F}_{\text{Tens}}$ are:

$ {\vec F}_{\text{Bend}} = \kappa\frac{\partial^4 {\vec r}}{\partial s^4}, $

(15)

$ {\vec F}_{\text{Tens}} = \beta\frac{\partial^2{\vec r}}{\partial s^2}. $

(16)

Applying this result to elastic light scattering, the model yields correct results in the flexible coil limit ^[108], but it fails at very high stiffness, deviating strongly from the theoretical solution obtained for rigid rods ^[109].

For computer simulations, a discrete polymer chain model has to be employed. Here, the bending rigidity, i.e., the stiffness of the chains which are composed of $N$ mass points (monomers), is incorporated into the CG polymer model by the following bending potential

$ U_{\text{bend}} = \frac{\kappa}{2}\sum\limits_{i=1}^{N-1} ({\vec u}_{i+1}-{\vec u}_i)^2 = \kappa\sum\limits_{i=1}^{N-1} (1-{\vec u}_{i+1}\cdot{\vec u}_i) = \kappa \sum\limits_{i=1}^{N-1} (1 - \cos\theta_{i, i+1})\;, $

(17)

where $\vec{u}$ is the unit bond vector ${\vec u}_i = {\vec r}_{i+1} -{\vec r}_i/{\mid {\vec r}_{i+1} -{\vec r}_i\mid }$, connecting consecutive monomers and $\vec{r}_i$ is the position vector to the $i$-th monomer. Parameter $\kappa$ scales the strength of the bending energy in the model. In the discrete numerical model, $L_p$ is fixed by parameter $\kappa$. A straightforward definition of the persistence length $L_p$ in a discrete simulation model with respect to the bond length $d_0$ is given by:

$ L_p = \frac{\kappa d_0}{k_BT}\;. $

(18)

In the simulations, $L_p$ is varied in the range of $L_p\in \{5\, d_0, 10\, d_0, 20\, d_0, 50\, d_0\}$, and is thus large compared to one bond length $d_0$, which is the smallest length scale on which the model chains still exhibit flexibility. This range of $L_p$ is also larger than the persistence lengths used in most previous MD studies of semiflexible chains which were in the range of $L_p\in \{1\; d_0, ..., 5\, d_0\}$ ^{[102,103,104,105]}. Monodisperse chains with monomer numbers $N=350$ and $N=700$ are simulated. Thus, the semiflexible regime of polymer chains for which

$ d_0 \ll L_p \ll L=(N-1)\, d_0 $

(19)

is investigated.

The simplest model for the description polymer chain dynamics in a melt, is the Rouse model ^[107]. It is a three parameter model, including the monomeric friction ($\xi$), the chain connectivity (modeled through harmonic springs with mean-squared bond length $d_0{}^2$) and the degree of polymerization of the chains $N$. This model is known to describe the dynamics of short, unentangled melts reasonably well, though deviations appear at monomeric length scales, which are affected by local excluded volume interactions and chain stiffness. For long, entangled chains, the Rouse model describes the dynamics at intermediate time/length scales even though the longer scale dynamics are strongly affected by constraints formed by surrounding chains. The Rouse modes, $p = 0, 1, 2, ..., N-1$, of a chain of length $N$ are defined as:

$ \vec{X}_p = \sqrt{\frac{2}{N}}\sum\limits_{i=1}^N \vec{r}_i\cos\left[\frac{p\pi}{N}\left(i-\frac{1}{2}\right)\right]. $

(20)

The $p=0$ mode describes the motion of the chain center-of-mass, while the modes with $p\leq 1$ describe internal relaxations in polymer chains with a mode number $p$ corresponding to a subchain of $(N-1)/p$ segments. Hence, a Rouse mode analysis is in essence analogous to a Fourier analysis.

Starting from the Langevin Equation of the Rouse model,

$ \xi\, \dot{\vec{r}}(s, t) = k\frac{{\partial}^2}{\partial s^2} \vec{r}(s, t) + {\vec{F}}^S(s, t), $

(21)

where $s$ again is some parameter along the polymer chain that measures length—for details, see ^[107]—the stiffness of the polymer chains can be taken into account by introducing an additional entropic bending force $\vec{F}^B$ of the form

$ \vec{F}^B(s, t) = -L_pk_BT\frac{\partial^4}{\partial n^4}\vec{r}(s, t), $

(22)

which can be derived by applying Hamilton's principle to Eq. (14) ^[110]. The equation of motion then becomes

$ \xi\frac{\partial}{\partial t}{\vec{r}}(s, t) = \frac{3k_BT}{2l_P}\frac{\partial^2}{\partial s^2}{\vec{r}}(s, t) - \frac{3k_BTL_p}{2}\frac{\partial^4}{\partial s^4}{\vec{r}}(s, t) + {\vec{F}}^S(s, t)\;, $

(23)

where $\vec{F}^S$ is a Gaussian stochastic force. Eq. (23) is similar to the equation of motion derived by Harris and Hearst ^[109] and by Harnau et al. ^[111] who, as in the Rouse model, assume a Gaussian distribution of each bond length. In essence, (23) is based on the Rouse model for flexible Gaussian chains and the introduced bending force can be considered as a small and local perturbation of the system. A solution of (23) can be achieved in terms of a normal Rouse mode analysis as shown in ^[112]:

$ \langle{\vec{\tilde X}}_p{}^2(0)\rangle = \frac{k_BT}{k_p{}^{\text{semi}}} = \left[\frac{3\pi^2}{LL_p}p^2 + \frac{3L_p\pi^4}{L^3}p^4 \right]^{-1}\;, $

(24)

where the effective force constant of semiflexible chains is introduced as:

$ k_p{}^{\text{semi}} = k_p{}^b + k_p = (3k_BTL_p\pi^4p^{4})/(L^3) + (3k_BT\pi^2p^2)/(LL_p)\;. $

(25)

The value $k_p{}^b $ is due to the mechanical bending force and $k_p$ arises from the entropic tension. The crossover scaling behavior of single polymers of different persistence lengths $L_p$ is displayed in Figure 17. As one can see, these chains follow the corresponding scaling laws as predicted from theory (Eq. (24)). Introducing a critical Rouse mode $p_c$ upon which the crossover scaling from a $p^{-2}$ to $p^{-4}$ behavior, i.e., from entropic tension modes to mechanical bending modes occurs, we can see, that all different simulations, independent of system size $N$ fall on one universal master curve as displayed in Figure 17. Since we discuss here CG models of polymers, much longer persistence lengths $L_p$ (see Eq. (18)) than in atomistic simulations of universal scaling properties can be considered. Atomistic simulations done by Paul et al. ^{[102,105,113]} claimed the scaling exponents to be in the range of $\approx p^{-2.3}$ to $\approx p^{-3.0}$, while theory predicts for the semiflexible regime a universal exponent of $p^{-4.0}$, where $p$ is the Rouse mode number of Eq. (24). In the simulations presented in ^[19] it was shown for the first time that the correct (theoretically predicted) crossover scaling in the transition from fully flexible to semiflexible behavior of linear macromolecules can be obtained.

---

5.5. Coarse-Graining of Lipid Bilayer Membranes

In this section, I discuss a recently introduced CG model of lipid molecules which self-aggregate to form stable bilayer membranes and which are subsequently exposed to shock waves. Biological cell membranes are very interesting and complex supramolecular structures which form a barrier between different compartments, or organelles, of cells. Membranes also exhibit many chemical reactions essential to the existence and functioning of a cell ^[114]. For example, the plasma membrane acts as a boundary that prevents the contents of a cell from escaping and mixing with the surrounding medium. Simultaneously, a cell's plasma membrane enables the passage of critical nutrients into the cell and the passage of waste products out. The most common membrane constituents are lipid molecules with two physically separated subdomains, usually an elongated hydrophobic domain made up of fatty acid tails, associated with a hydrophilic head group. The most abundant lipids in cell membranes are the phospholipids, molecules in which the hydrophilic head is linked to the rest of the lipid through a phosphate group ^[90], cf. Figure 18. The hydrophilic head groups of lipids dissolve readily in water because they contain polar groups that are easily incorporated in the hydrogen-bonding network of the surrounding water. In contrast, the hydrophobic hydrocarbon tails are uncharged and non-polar and thus try to aggregate in energetically and entropically favourable structures that minimize their contact with surrounding water molecules. The nearly cylindrical shape of most membrane lipids makes the bilayer the most common geometrical organization for spontaneous self-assembly of lipid molecules in aqueous environment.

All-atom MD simulations of lipid bilayers which resolve the dynamics of individual atoms are limited to fairly small membrane samples (tens of nanometers in extension) and very short time scales of at most a few hundred nanoseconds ^[116], leaving CG simulations as the only currently available computational tool to access mesoscale phenomena such as the dynamics of molecular self-assembly of lipid molecules.

There is a very large body of literature on computational studies of the static and dynamic properties of biomembranes using atomistic and CG modeling approaches, see e.g., ^[117,118] which have been reviewed in depth, e.g., by Pandit and Scott ^[119] and Woods and Mulholland ^[95].

With the rise of so-called solvent-free, or implicit simulation schemes for membranes, which became fashionable at the turn of the millennium, the number of publications in this field has constantly increased. Solvent-free models of lipid bilayer structures do not explicitly take the fluid molecules of the aqueous environment into account. These models are either based on a stochastic Langevin equation of motion accounting for the Brownian random motion of the fluid molecules, or on the complete modeling of all effects of the solvent by effective interaction potentials between the constituent particles of the membrane only ^{[120,121,122,123]}. As one is usually only interested in the structural and dynamic details of the membrane and not in the surrounding fluid, this allows to reduce the number of necessary integrations and thus the computational costs considerably. Many of the used potentials are either derived from intuition, from standard potentials routinely being using in polymer physics or from quantum chemical calculations of force field parameterizations ^{[117,124,125,126]}.

---

5.5.1. Phase Diagram of a Lipid Bilayer Membrane Model

When employing a solvent-free model, the hydrophobic interactions of the lipid tails with water can be considered using an attractive potential between all tail particles of different lipid chains as indicated in Figure 18:

$ $$\begin{equation} \label{eq:pot4} \Phi_{\text{attr}}(r) = \left\{ \begin{array}{ccc} -\xi &\quad\mathrm{for}& r<r_{\text{cut}}\, , \\ -\xi\cos^2\left(\frac{\pi\left(r-r_c\right)}{2h_c}\right) &\quad\mathrm{for}& r_c\leq r\leq r_c + h_c\, , \\ 0 &\quad\mathrm{else}& \, , \end{array} \right. \end{equation}$$ $

(26)

where $h_c$ determines the range of the attractive potential, i.e., the larger $h_c$ the larger the hydrophobicity of the lipid tails, and $r_{\text{cut}}=\sqrt[6]{2}$. A tunable potential of this type (involving a cosine function) was introduced by Steinhauser in polymer physics for the realistic MD simulation of polymer-solvent interactions with varying solvent qualities ^[83]. The decay range $h_c$ of Eq. (26) directly determines the bending rigidity and rupture strength of the membrane. Thus, a variation of this parameter in a computational study can be used to explore the properties of self-assembled lipid bilayers with a broad range of mechanical properties. One finds stable bilayer structures within the interval $h_c=\left[0.8 \: \sigma, 1.4 \: \sigma\right]$ at a temperature $T=310 \: K$, which corresponds to physiological conditions. For $h_c \lesssim 0.7 \: \sigma$, the bilayer dissolves, while it is crystalline for $h_c \gtrsim 1.5 \: \sigma$.

In Figure 19, the phase diagram of lipid membranes based on the two parameters temperature $T$ and attractive potential range $h_c$ is displayed. Each data point was actually generated from averaging and analyzing 20 different simulations for the specific combination of $(T, h_c)$. The fluid character of the bilayer within the stable range is verified by computing the in-plane diffusion coefficient $D=0.03 -0.06 \: \sigma^2 / \tau$ and observing a non-zero flip-flop rate, i.e., the probability per unit time that a single lipid changes from one surface to the other. The phase diagram of membrane structures obtained with this model was calculated by averaging $100$ simulation runs for each displayed combination of temperature $T$ and attractive potential range $h_c$.

---

5.6. Simulation of Shock Wave Damage in Coarse-Grained Models of Membranes

CG models are nowadays routinely used in polymer and biophysics, mostly for equilibrium ^{[127,128,129,130,131,132,133]}, but also for non-equilibrium studies ^{[75,134,135,136]} of lipid bilayer membranes and provide a description of reduced complexity with respect to the molecular degrees of freedom ^[72,137,138].

Practically all reported simulation studies of the static and dynamic properties of biomembranes have been performed very close to equilibrium. However, the exploration of the interaction of shock waves with biomembrane structures constitutes a non-equilibrium physical process. In the few existing computational studies exploring shock waves in soft matter none takes properly into account the particular physical conditions of a shock wave and they are limited to unrealistically small systems with explicitly modeled solvent ^{[139,140,141]}.

In a recent research paper, I suggested a new multiscale approach to this problem, discussing the fact that even in the largest existing all-atom simulation published even to date ^[139,140], the size of the considered system was several orders of magnitude too small in size and too limited with regards to the time scale to be able to capture any relevant effects of shock waves in membranes observed in experiments, such as transient permeability and subsequent self-repair of parts of the membrane in an eukariotic cell ^[134]. It is well known that the exposure of biological cells to shock waves can cause damage of varying extend to the cell membrane ^{[142,143,144,145]} depending on the pressure level of the shock wave ^[5,91]. However, the exact mechanisms by which shock waves interact with membranes of biological cells is vastly unknown.

In medicine, the use of shock waves as a form of therapy has two major areas of application:

1. Gene therapy: Here, the plasma membrane is transiently punctured by the treatment with shock waves, which can be used for delivering genetic material or therapeutic molecules into the cytoplasm ^[8].

2. High Intensity Focused Ultrasound (HIFU): Here, tumor tissue is treated with HIFU which transports energy to a focal point where it is converted into heat, leading to denaturation of proteins in the cells (coagulative necrosis). A secondary and hard to control destructive effect during HIFU treatment is caused by acoustic cavitation, i.e., the formation and implosion of small gas bubbles at the focal point.

Experimentally, direct visualization of the dynamics of membrane rupture is very difficult to achieve, because a typical cell membrane has thickness $d \approx 5 \: \mathrm{nm}$, while the pressure front of a shock wave travels at supersonic speeds, i.e., with a velocity $v_s \gg 1430 \: \mathrm{ms^{-1}}$ in water. The time-scale during which a shock wave interacts with a cell membrane is thus on the order of a few picoseconds.

The intrinsic structural complexity of cell membranes bears an additional computational challenge: A membrane is essentially a thin fluid layer with immersed proteins and carbohydrates, stabilized by a subtle balance of competing hydrophobic and hydrophilic interactions which have to be reflected in a corresponding CG model.

---

5.6.1. Membrane Damage: Effects of Shock Wave Speed

Using a DPDE thermostat ^[146], recent shock wave computer experiments on square lipid bilayer patches were performed for values of the lipid interaction parameter $h_c$ of Eq. (26) in the interval $h_c \in \left[0.8, 1.4\right]\sigma$, see Figure 20. For $h_c \lesssim 0.7 \: \sigma$, the bilayer is just instable and dissolves, while it attains a crystalline shape for $h_c \gtrsim 1.5 \: \sigma$. Hence, the investigated range of interactions is representative of a wide spectrum of lipid bilayers with different mechanical stability. Initial configurations are taken from equilibrium runs performed in a tensionless state. The here shown system consists of $N_{\text{lip}}=3872$ lipids, spread on an area of $A =15.01^2 \: \mathrm{nm}^2$. The simulation box size along the direction normal to the bilayer is chosen as $L=63.08 \: \mathrm{nm}$ The total number of particles is $N \approx 1.1\times 10^6$. The driving speed of the piston that generates the shock wave is varied in the range $v_p \in \left[1892,5676\right]\mathrm{ms}^{-1}$ in order to produce shock waves with different supersonic velocities. After 200 time steps of $\delta t=2.87 \: \mathrm{fs}$, the piston is stopped, while the initiated shock wave continues to travel further along $L$. The lipid bilayer is placed $9.435 \: \mathrm{nm}$ in front of the initial piston position, far enough away not to be hit by the piston. Snapshots of an exemplary simulation are shown in Figure 21, where the shock wave was initiated by a piston with velocity $v_p=4730 \: \mathrm{m}\mathrm{s}^{-1}$.

There are a number of different phenomena associated with lipid bilayer damage:

• The bilayer can tear along well localized paths leaving the majority of its area intact;

• The two bilayers can inter-penetrate upon compression;

• Single lipids can move out of their equilibrium positions and orientations via diffusive mechanisms.

All of these effects may conveniently be combined into a single scalar order parameter, $\Psi$, based on the projection of the particle pair distribution function on rotational invariants ^[148]

$ C(r) = \frac{1}{ \pi N}\, \sum\limits_{i > j}^N \delta(r - r_{ij})\, \frac{\vec{e}_i \cdot \vec{e}_j}{r^2}, $

(27)

where the sum runs over all pairs of lipids $i$ and $j$, $r_{ij}$ is the pair distance as measured between the central beads, and $\vec{e}_i$ and $\vec{e}_j$ are orientation vectors (the normalized distance vectors between the first and the third bead of each lipid). $C(r)$ is similar to the familiar radial distribution function, except that it is additionally weighted by relative orientations. To convert this distance dependent correlation function into a single scalar, $C(r)$ is integrated up to a distance $r_c=1.0 \: \mathrm{nm}$, which is chosen such as to include the first coordination shell, i.e., the first peak in $C(r)$, in the undamaged bilayer structure:

$ \Psi = 4 \pi \int\limits_{0}^{r_c}\, \mathrm{d}r\, C(r) r^2\, . $

(28)

Disruption of the equilibrium bilayer configuration affects this order parameter in two different ways: If the mutual orientation of neighboring lipids deviates from parallel alignment, or if lipids are separated from each other beyond their equilibrium distance $r_c$, the value of $\Psi$ is reduced, thus providing a quantitative means to characterize lipid bilayer damage.

Figure 22 shows the damage caused by shock waves of different impact velocities $v_s$ for membrane systems of medium stability with an interaction parameter $h_c=1 \: \sigma$. As one can see, higher impact velocities cause a stronger decrease of $\Psi$, i.e., larger damage. The initial descent of $\Psi$ is almost instantaneous for the here simulated shock front speeds, however, for large simulation times $\Psi$ increases within a few ps, indicating reversible damage. This is consistent with the experimentally observed self-repair of biomembranes that are transiently permeable when exposed to shock waves ^[4]. For small shock front velocities $v_s \lesssim 3050 \: \mathrm{m}\mathrm{s}^{-1}$, no such clear increase of $\Psi$ is observed. Thus, the existence of a critical shock front velocity for self-recovery of membranes can be deduced from our calculations.

---

6. Conclusions

The combination of shock wave physics, coarse-graining and multiscale modeling in computer simulations will continue to be a very important part of state-of-the-art computational materials science. The combination of research results and expertise of these different research areas has led to the new emerging area of investigating shock wave effects in biological, i.e., soft matter systems and will eventually overcome (at least partly) the traditional and artificial barriers between scientists from different branches such as mathematics, physics, computer science, biology, chemistry, etc. who usually do not work well together. The importance of materials modeling across length and time scales will likely further increase in the future with the development of novel computational techniques and the advent of exascale computers. Currently, mostly sequential multiscale modeling and bridging of few scales is state of the art, while concurrent modeling concepts begin to emerge. It can be expected that the future will bring more concurrent modeling for a variety of materials that would enable design of materials at the macroscale by using an inverse engineering design methodology. Concurrent modeling can be achieved through further development of a systematic coarse-graining methodology that retains microstructure heterogeneities and chemistry of the system at higher scales while the reverse process will utilize general backmapping algorithms that can uncover underlying atomistic structure and local heterogeneity from the high-level coarse-grained representation.

New algorithms will likely improve the efficiency and validity of computations at every scale by making the representation of the total material system more accurate with better estimates of the uncertainties that have arisen. For example, embedded and adaptive methods may enable more accurate calculations of the material system at the various scales (quantum, atomistic, micro-and meso-scales) subjected to a dynamically adjusted environment. Current, predetermined a priori representation of the multiscale material system for atomistic simulations (force fields) or for FE techniques (constitutive equations) could be calculated on the fly, thus adding accuracy and less uncertainty to the material system simulation.

Undoubtedly, with increased coupling of different algorithms across spatiotemporal scales, understanding deficiencies and the range of applicability of models, as well as estimation of error on simulations, will become of vital importance. The field of verification, validation, and uncertainty quantification will become an integral part of concurrent multiscale modeling. In the future, we can expect an increasing integration of multiphysics within multiscale modeling, allowing for realistic simulations under various external fields and extreme conditions. Bulk calculations might be replaced with realistic simulations including local heterogeneities, such as stochastically distributed defects, grain and interface boundaries or microphase-separated morphologies that can dynamically evolve upon influence of external fields.

What would probably really be needed to address the interesting developments in computational science outlined in this review is a change of traditional curricula at universities worldwide that barely offer inter-or multi-disciplinary courses that combine e.g., physical theories and mathematical modeling with traditional wet-lab cancer research or e.g., look into the application of quantum mechanics in biology—another new emerging area that has become known under the name quantum biology. Also, computational science courses that combine techniques and research results from different areas in integrated courses are currently rarely offered. One problem here might be, that the present academic teachers at universities themselves have not grown up in a multi-and interdisciplinary environment when they obtained their own degrees, and then most likely have been busy for many years following the infamous "publish-or-perish" fetish which most likely produces countless numbers of redundant research papers that serve mostly only for one thing: To act as research currency for academic careers and funding.

---

Acknowledgments

The author acknowledges financial support from the German Federal Ministry for Economic Affairs and Energy under grant no. 50 LZ 1502.

---

Conflict of Interest

The author declares that there are no conflicts of interest related to this review article.

---