Parking-garage structures in astrophysics and biophysics
C. J. Horowitz,
1,
∗
D. K. Berry,
1,
†
M. E. Caplan,
1,
‡
Greg Huber,
2,
§
and A. S. Schneider
3,
¶
1
Center for Exploration of Energy and Matter and Department of Physics,
Indiana University, Bloomington, IN 47405, USA
2
Kavli Institute for Theoretical Physics, Kohn Hall,
University of California, Santa Barbara, California 93106-4030, USA
3
TAPIR, California Institute of Technology, Pasadena, CA 91125, USA
(Dated: September 2, 2015)
A striking shape was recently observed for the cellular organelle endoplasmic reticulum consisting
of stacked sheets connected by helical ramps [1]. This shape is interesting both for its biological
function, to synthesize proteins using an increased surface area for ribosome factories, and its geo-
metric properties that may be insensitive to details of the microscopic interactions. In the present
work, we find very similar shapes in our molecular dynamics simulations of the nuclear pasta phases
of dense nuclear matter that are expected deep in the crust of neutron stars. There are dramatic
differences between nuclear pasta and terrestrial cell biology. Nuclear pasta is 14 orders of magni-
tude denser than the aqueous environs of the cell nucleus and involves strong interactions between
protons and neutrons, while cellular scale biology is dominated by the entropy of water and com-
plex assemblies of biomolecules. Nonetheless the very similar geometry suggests both systems may
have similar coarse-grained dynamics and that the shapes are indeed determined by geometrical
considerations, independent of microscopic details. Many of our simulations self-assemble into flat
sheets connected by helical ramps. These ramps may impact the thermal and electrical conductiv-
ities, viscosity, shear modulus, and breaking strain of neutron star crust. The interaction we use,
with Coulomb frustration, may provide a simple model system that reproduces many biologically
important shapes.
PACS numbers: 26.60.-c,02.70.Ns, 87.10.Tf, 87.15.ap, 87.16.D-,87.16.Tb
Membrane-bound cellular organelles have characteris-
tic shapes, with the endoplasmic reticulum (ER) being
particularly striking. The ER is an extensive organelle
displaying three distinct, yet connected, morphologies:
tubes, sheets and the spherical envelope around the cell
nucleus. Recent advances in serial sectioning and electron
microscopy have revealed the stacked ER sheets to be
connected by helical structures [1]. Just as spiral ramps
connect the levels of a multilevel parking garage, these
so-called “Terasaki ramps” allow sheets to connect yet re-
main parallel over scales large relative to the membrane
thickness. This parking-garage shape is interesting for
both its biological function, to synthesize proteins using
an increased surface area for ribosome factories, and its
mathematical property as the minimizer of a geometrical
Hamiltonian largely insensitive to details of the micro-
scopic interactions [2]. In these models of membrane me-
chanics, the Terasaki ramps are stable, topological struc-
tures, akin to screw dislocations, that affect the entire
morphology of an organelle critical to the metabolism of
eukaryotic cells.
Independent of considerations from terrestrial biology,
in this Letter we find very similar shapes arising in our
∗
Electronic address: horowit@indiana.edu
†
Electronic address: dkberry@iu.edu
‡
Electronic address: mecaplan@indiana.edu
§
Electronic address: huber@kitp.ucsb.edu
¶
Electronic address: andschn@indiana.edu
molecular dynamics (MD) simulations of dense nuclear
phases. These phases are expected in neutron-star crusts
and in core-collapse supernovae [3]. If matter is com-
pressed to near nuclear densities, competition between
short-range nuclear attraction and long-range Coulomb
repulsion can lead to a variety of complex shapes, so-
called “nuclear pasta” [4, 5]. Note that Groenewold
et
al.
[6] discuss colloidal cluster phases arising from short
range attraction and coulomb repulsion. For a particu-
lar density and proton fraction, we also find flat sheets
connected by spiral ramps. In a neutron star, electrons
scattering from the spiral ramps reduce both the thermal
and electrical conductivity. This may impact X-ray ob-
servations of crust cooling in transiently accreting stars
[3], and may also lead to the decay of neutron-star mag-
netic fields after about a million years [7].
There are, to be sure, dramatic differences between
nuclear pasta and terrestrial cell biology. Nuclear pasta
has a density near 10
14
g/cm
3
, fully fourteen orders of
magnitude more dense than the aqueous environs of the
cell nucleus. Furthermore, nuclear pasta involves strong
interactions between neutrons and protons in addition
to electromagnetic interactions, while cellular-scale bi-
ology is highly-screened, highly-overdamped and domi-
nated by the entropy of water and complex assemblies
of biomolecules. Nonetheless the strikingly similar ge-
ometry suggests both systems may have similar coarse-
grained dynamics.
In nuclear physics, the semi empirical mass formula
predicts the binding energy
BE
of a nucleus with
A
nu-
arXiv:1509.00410v1 [nucl-th] 30 Aug 2015
2
cleons (neutrons plus protons) and
Z
protons [8, 9],
BE
=
a
v
A
−
a
s
A
2
/
3
−
a
c
Z
2
/A
1
/
3
...
(1)
Here
a
v
,
a
s
,
a
c
are constants describing volume, surface,
and Coulomb energies. In addition there are other con-
tributions from the symmetry energy and pairing that
will not be important here. Competition between surface
and Coulomb energy contributions can lead to complex
shapes. We perform simulations at a density of
n
= 0
.
05
fm
−
3
. This corresponds to a packing fraction of 5/16 of
nuclear saturation density,
n
0
= 0
.
16 nucleons per fm
3
.
Here the system forms flat sheets (lasagne) that are con-
siderably thicker than the size of a single nucleon. This
thickness is determined from a balance of surface and
Coulomb energies. Note that there are a variety of dif-
ferent shapes with almost the same energy. For example,
if the density is decreased somewhat, the system forms
rods (spaghetti) instead of flat sheets and at still lower
densities the system forms spheres representing isolated
nuclei.
Fluctuations about these simple shapes could have
low excitation energies that may depend on sub-domaint
terms in Eq. 1. Reinhard
et al.
[12], see also [13], use a
leptodermous (thin-skinned) expansion and density func-
tionals in order to calculate a curvature energy term for
Eq. 1 that goes like
a
curve
A
1
/
3
. The impact of such an
A
1
/
3
term on nuclear pasta shapes was considered in ref.
[14]. The curvature energy could be important for nu-
clear fission where a nucleus dramatically changes shape.
However, Reinhard
et al.
find they need to calculate the
energy of very large nuclei with
A
in the thousands in
order to extract
a
curve
theoretically. Therefore, this co-
efficient may not be reliably determined from measured
nuclear binding energies that are known only for a limited
range of
A
.
Instead, we discuss a very different approach that
has been employed previously in biophysics, but not
(to our knowledege) in nuclear physics. We consider a
Helfrich-Canham Hamiltonian
H
0
[10, 11] that involves
a quadratic form in the surface principal curvatures
C
1
and
C
2
including both the mean curvature (
C
1
+
C
2
)
/
2
and Gaussian curvature
C
1
C
2
,
H
0
=
1
2
B
∫
dS
(
C
1
+
C
2
)
2
+
̄
B
∫
dSC
1
C
2
.
(2)
Here
∫
dS
is an integral over the surface area and
B
and
̄
B
are constants representing effective rigidity mod-
uli. This energy functional was applied in refs. [1, 2] to
biological systems. Our justification for also discussing
nuclear pasta shapes with the help of Eq. 2 is determined
only after the fact. This equation predicts spiral ramp
shapes, and their arrangements, very similar to what we
find in our MD simulations, see below.
One relatively low energy solution for Eq. 2 is lasagne
with flat surfaces where
C
1
=
C
2
= 0. This gives
H
0
= 0
.
Another solution involves spiral ramps with
C
1
=
−
C
2
.
Here the mean curvature is still zero, but the Gaussian
curvature is negative so that
H
0
<
0. Thus the Gaussian
curvature term may stabilize spiral ramp configurations
in both biological membranes and nuclear pasta. This
term may also stabilize configurations that have addi-
tional holes such as the “nuclear waffle” shapes found
in ref. [15]. These shapes consist of flat sheets with a
two dimensional array of holes. Similarly, the ER dis-
plays a waffle-type pattern which has been referred to
in the biological literature as “fenestrations” [16]. Al-
ternatively there may be torus or donut-shaped super-
heavy nuclei [17] where the donut shape both reduces
the large Coulomb energy and is further stabilized by
the Gaussian-curvature term.
In this paper we study the self-assembly of these spi-
ral ramp configurations with molecular dynamics simula-
tions of a simple (semi)classical model of nuclear matter.
Our simulations are for nuclear pasta but they may also
have implications for phospholipid bilayer membranes.
Our MD formalism is the same as that used by
Horowitz
et al.
in previous works [3, 18–22] and is briefly
reviewed here. It is very similar to a model used by others
[23, 24]. Our simulation volume is a cubic box with peri-
odic boundary conditions which contains point-like pro-
tons and neutrons with mass
M
= 939 MeV. Electrons
are assumed to form a degenerate relativistic Fermi gas
and are not explicitly included in the simulations. Pro-
tons and neutrons interact via the two-body potentials:
V
np
(
r
) =
a
e
−
r
2
/
Λ
+
b
e
−
r
2
/
2Λ
(3a)
V
nn
(
r
) =
a
e
−
r
2
/
Λ
+
c
e
−
r
2
/
2Λ
(3b)
V
pp
(
r
) =
a
e
−
r
2
/
Λ
+
c
e
−
r
2
/
2Λ
+
α
r
e
−
r/λ
.
(3c)
The
n
and
p
indices indicate whether the potential de-
scribes a neutron-proton, a neutron-neutron, or a proton-
proton interaction. Meanwhile,
r
is the separation be-
tween each pair of interacting nucleons and quantities
a
= 110 MeV,
b
=
−
50 MeV,
c
=
−
2 MeV, and Λ = 1
.
25
fm
2
are parameters of the model. Their values were cho-
sen in ref. [18] to approximately reproduce some bulk
properties of pure neutron matter and symmetric nuclear
matter, as well as the binding energies of selected nuclei.
The screening length
λ
is chosen to be 10 fm. Equation
3 describes an intermediate range attraction between
p
and
n
which binds nuclei and then a short range repul-
sion that causes nuclear matter to saturate at a density
n
0
= 0
.
16 fm
−
3
. Finally there is a long range Coulomb
repulsion between protons. All of the simulations in this
paper are at a density of
n
= 0
.
05 fm
−
3
, a composition of
40% protons, 60% neutrons, a fixed temperature
kT
= 1
MeV, and use a MD time step of 2 fm/c.
To study self-assembly of these ramps, we start from
initial conditions where the particle positions are uni-
formly distributed in the simulation volume with a Boltz-
mann velocity distribution. Figure 1 shows the configura-
tion of a 40,000 particle simulation at times from 40,000
to 2,000,000 fm/c. The proton density
n
p
is shown where
3
FIG. 1: (Color online) Self-assembly of a parking garage struc-
ture in a MD simulation with 40,000 nucleons that started
from uniform random initial positions. Shown in panels (a)
through (f) are the configuration after simulation times of
40,000, 400,000, 800,000, 1,200,000, 1,600,000 and 2,000,000
fm/c respectively. The color map at right is discussed in the
text. Figure produced with Paraview [25].
each proton is assumed to have a Gaussian shape of 1
fm
3
volume [15] and the opacity of the colorscale is 0 for
n
p
= 0
.
00 to 0.02 fm
−
3
and then increases linearly to
1 at
n
p
= 0
.
04 fm
−
3
. A light cream color corresponds
to high density sheets, while lower density surfaces are
shown in brown. This system undergoes the following
self-assembly steps: (1) the low density system collapses
locally to form higher density filaments that meet in junc-
tions, see Fig. 1 (a). This includes the formation of a
number of topological holes. (2) Next, the filaments start
to grow to form curved sheets, Fig. 1 (b-c), (3) these
sheets then start to straighten out over longer length
scales, Fig. 1 (d). (4) Boundaries between “domains”
of sheets with different orientation form four left-handed
and four right handed helical ramps, see Table I. The
sheets straighten out over the full simulation volume, see
Fig. 1 (e-f), and the ramps move together to form the
dipole pattern shown in Fig. 2 (a). This pattern has
four left-handed ramps to the left and four right-handed
ramps to the right. For this dipole pattern, the ramps
are seen in Fig. 2 (c) to make about a 45 degree angle
with the flat sheets. We find this final configuration to
be stable for times of at least 10,000,000 fm/c.
To study the dependence on boundary conditions, we
perform simulations with different numbers of particles
and correspondingly different sized simulation volumes,
see Table I. The smallest simulation, with only 20,000
particles, forms uniform flat sheets without any spiral
ramps (not shown). A 50,000 particle simulation forms
four left handed and four right handed ramps as shown
in Fig. 2 (e). This is also a dipole pattern with the left-
handed ramps to the left and the right-handed ramps to
the right. As a result the ramps make a 45 degree angle
with the sheets, see Fig. 2 (f). Finally a 75,000 parti-
cle simulation forms four ramps in a quadrupole pattern
where one left-handed and one right-handed ramp are to
FIG. 2: (Color online) Three axis views, X (left), Y(center),
and Z (right), of the final configuration of three MD simula-
tions. Panels (a)-(c) at top are for the 40,000 nucleon sim-
ulation shown in Fig. 1. Panels (d)-(f) at center are for a
simulation with 50,000 nucleons and panels (g)-(i) at bottom
are for a simulation with 75,000 nucleons. The handedness of
the helical holes in panels (a), (e), and (h) are indicated with
L for left-handed and R for right-handed.
the left as shown in Fig. 2 (h). For a quadruple pattern
the ramps are observed to make a 90 degree angle with
the sheets, see Fig. 2 (i) and Table I.
The arrangement of the helical ramps shown in Fig.
2 panels (a), (e) and (h) agrees very well with theoret-
ical predictions in refs. [1, 2] that are based on Eq. 2.
Guven
et al.
argue that tension in the sheets leads to
an effective long range attraction between two ramps of
opposite chirality that draws them together until a short
range repulsive bending force stabilizes the pair of ramps
at a characteristic distance. Some features of our simula-
tions that are in agreement with Guven
et al.
are: (1) we
find an equal number of left-handed (L) and right-handed
(R) ramps, (2) the ramps are all relatively close together
with a similar characteristic spacing, and (3) ramps in a
quadrupole pattern make a 90 degree angle with respect
to the sheets while ramps in a dipole pattern make an
approximately 45 degree angle with the sheets. These
common features suggest, after the fact, that Eq. 2 may
also apply to our nuclear pasta model system.
Most of our simulations form flat sheets connected by
helical ramps. However we are able to obtain only flat
sheets, without any ramps, if we add a small one-body
potential to the system for early times that biases the
formation of only flat sheets, or start a simulation in a
smaller box at high densities where the system is nearly
4
TABLE I: Number of ramps
N
r
in MD simulations with
different numbers of particles
N
.
NN
r
Pattern of ramps
Angle of ramps (deg)
20,000
0
40,000
8
Dipole
45
50,000
8
Dipole
45
75,000
4
Quadrupole
90
uniform and then very slowly expand the box during the
simulation until the system reaches the same final density
of
n
= 0
.
05 fm
−
3
.
FIG. 3: (A) Scanning electron micrograph of a thin slice of
sheet-like ER from a mouse salivary gland. Serial section-
ing allows three-dimensional reconstruction from large num-
bers of cross-sections (scale bar = 200 nm). (B) 3D recon-
struction of a left-handed Terasaki ramp that appears in the
black-outlined region in (A). (C) 3D reconstruction of a right-
handed Terasaki ramp that appears in the white-outlined re-
gion in (A). (Figure adapted from Terasaki et al. [1])
We compare our results to biological observations.
Reconstructing the three-dimensional geometry of ER
sheets from two-dimensional serial sections reveal that
the continuity of parallel sheets comes about through
the helical winding of the “exposed” sheet edge through
space,
i.e.
the Terasaki ramp, see Fig. 3 . The core of
the ramp (cytosolic side of the membrane) has a highly
negative Gaussian curvature, but potentially small mean
curvature, given the opposite signs of the two principal
radii of curvature. No considerations seem to set a pre-
ferred handedness to the ramps, and the observations
at hand, though statistically small, are consistent with
right- and left-handed ramps in equal numbers through-
out the organelle. If the sheet edges are treated as ef-
fectively one-dimensional defects (thought to be stabi-
lized by membrane proteins on the cytosolic side), then
it is sufficient to treat the rest of the membrane by the
Helfrich-Canham Hamiltonian. A class of solutions min-
imizing that functional are minimal surfaces, which have
zero mean curvature and locally minimize area. One con-
sequence of the analysis in [2] is that whereas a single
Terasaki ramp has a logarithmically diverging energy, a
left-right dipole pair has a finite energy, and a double
pair of left-right-left-right ramps (a quadrupole) mini-
mizes it further. The existence of ramp dipoles is natural
in the biological context, though convincing evidence for
tightly correlated left-right pairs, or for pairs of pairs, is
currently lacking. As previously discussed, we find both
left-right pairs and pairs of pairs in our MD simulations
of nuclear pasta.
In conclusion, we have performed molecular dynam-
ics simulations using a simple classical model of nuclear
pasta. Many of our systems spontaneously self-assemble
to form flat sheets connected by helical ramps. This
geometry is very similar to that observed in the three-
dimensional structure of the sheet-like endoplasmic retic-
ulum. Seeing the same helical shapes in the extraordi-
narily different systems of nuclear pasta in neutron stars,
and in the membranes of eukaryotic cells strongly sug-
gests that these shapes follow from common geometric
considerations and may be independent of details of the
microscopic interactions.
The curvature energy could help stabilize these heli-
cal ramps, suggesting that they may be common in nu-
clear pasta. If true, electron scattering from the ramps
could reduce the electrical and thermal conductivity of
neutron-star crust. This may impact crust cooling [3]
and possibly lead to the decay of magnetic fields [7]. Al-
ternatively the helical ramps, that connect flat sheets,
may allow paired protons to percolate throughout the
system and make it superconducting. In addition, the
ramps may restrict how the sheets can move past each
other. This could lead to complex viscoelastic behavior
for nuclear pasta that could dampen r-mode oscillations
in rapidly rotating neutron stars, see for example [26].
Furthermore, the helical ramps likely will increase the
shear modulus of neutron-star crust and the frequencies
of crust oscillation modes. These modes may have been
observed as quasi periodic oscillations during magnetar
giant flares [27]. Finally, the ramps likely increase the
breaking strain, or mechanical strength, of neutron star
crust. As a result the crust can support larger moun-
tains that, on a rapidly rotating neutron star, efficiently
radiate gravitational waves [28].
One can study soft condensed-matter analogs of nu-
clear pasta in the laboratory even if direct experiments
on nuclear pasta are not feasible. This could provide
unique insights. For example, it may be very difficult
to predict the actual density of ramps in nuclear pasta
using only first principle simulations. Instead one can ob-
serve the actual density and pattern of Terasaki ramps
in analog laboratory systems.
Our simple interaction has Coulomb frustration with
short-ranged attraction and long-ranged repulsion. This
seems to provide a simple model system that reproduces
many biologically important shapes. One reason is the
fluid nature of bilayer membranes: interactions between
phospholipid molecules are determined to a large extent
by the repulsive term.
Computational advances have made and will make very
large-scale MD simulations for such systems “easy” and
these simulations will likely exhibit very rich varieties of
shapes and phases. Uncovering similarities between dis-
parate physical systems allows connections to be made
at the deeper level of symmetry, excitations and the ge-
ometry of topological defects.
5
Acknowledgments
We thank Gerardo Ortiz and Sima Setayeshgar for
helpful comments.
Part of this work was completed
at the Aspen Center for Physics, which is supported
by National Science Foundation grant PHY-1066293.
This research was supported in part by DOE grants
DE-FG02-87ER40365 (Indiana University) and DE-
SC0008808 (NUCLEI SciDAC Collaboration), and NSF
grant PHY11-25915 (KITP, UCSB).
[1] M. Terasaki
et al.
, Cell
154
, 285 (2013).
[2] J. Guven, G. Huber, and D. M. Valencia, Phys. Rev.
Lett.
113
, 188101 (2014).
[3] C. J. Horowitz, D. K. Berry, C. M. Briggs, M. E. Caplan,
A. Cumming, and A. S. Schneider, Phys. Rev. Lett.
114
,
031102 (2015).
[4] D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys.
Rev. Lett.
50
, 2066 (1983).
[5] M. Hashimoto, H. Seki, and M. Yamada, Progress of The-
oretical Physics
71
, 320 (1984).
[6] J. Groenewold, W. K. Kegel, J. Phys: Condens. Matter
16
, S4877 (2004).
[7] J. A. Pons, D. Vigano, and N. Rea, Nat. Phys.
9
, 431
(2013).
[8] C. F. von Weizacker, Zeitschrift fur Physik
96
, 431
(1935).
[9] P. Moller, W. D. Myers, W. J. Swiatecki, and J. Treiner,
Atomic data and nuclear data tables,
39
, 225 (1988).
[10] W. Helfrich, Z. Naturforsch. C
28
, 693 (1973).
[11] R. Canham, J. Theor. Biol.
26
, 61 (1970).
[12] P.-G. Reinhard, M. Bender, W. Nazarewicz, and T.
Vertse, Phys. Rev. C
73
, 014309 (2006).
[13] M. Durand, P. Schuck, X. Vifias, Z. Phys. A
346
, 87
(1993).
[14] K. Nakazato, K. Iida, and K. Oyamatsu, Phys. Rev. C
83
, 065811 (2011).
[15] A. S. Schneider, D. K. Berry, C. M. Briggs, M. E. Caplan,
and C. J. Horowitz, Phys. Rev. C
90
, 055805 (2014).
[16] M. Puhka
et al.
, Mol. Biol. Cell
23
, 2424 (2012).
[17] W. Nazarewicz, M. Bender, S. Cwiok, P.H. Heenen, A.T.
Kruppa, P.-G. Reinhard, T. Vertse, Nuclear Physics A
701
, 165c (2002).
[18] C. J. Horowitz, M. A. Perez-Garcia, and J.Piekarewicz,
Phys. Rev. C
69
, 045804 (2004).
[19] C. J. Horowitz, M. A. Perez-Garcia, J. Carriere, D. K.
Berry, and J. Piekarewicz, Phys. Rev. C
70
, 065806
(2004).
[20] C. J. Horowitz, M. A. Perez-Garcia, D. K. Berry, and J.
Piekarewicz, Phys. Rev. C
72
, 035801 (2005).
[21] C. J. Horowitz and D. K. Berry, Phys. Rev. C
78
, 035806
(2008).
[22] A. S. Schneider, C. J. Horowitz, J. Hughto, and D. K.
Berry, Phys. Rev. C
88
, 065807 (2013).
[23] C. O. Dorso, P. A. Gimenez Molinelli, and J. A. Lopez,
Phys. Rev. C
86
, 055805 (2012).
[24] P. N. Alcain, P. A. Gimenez Molinelli, J. I. Nichols, and
C. O. Dorso, Phys. Rev. C
89
, 055801 (2014).
[25] A. H. Squillacote, in The ParaView Guide: A Parallel Vi-
sualization Application (Kitware Inc., Clifton Park, New
York, 2007), p. 366.
[26] N. Andersson, K.D. Kokkotas, Int. J. Mod. Phys. D
10
,
381 (2001).
[27] A. L. Watts, T. E. Strohmayer, Advances in Space Re-
search,
40
, 1446 (2007).
[28] C. J. Horowitz and Kai Kadau, Phys. Rev. Lett.
102
,
191102 (2009).