of 5
PHYSICAL REVIEW C
94
, 055801 (2016)
“Parking-garage” structures in nuclear astrophysics and cellular biophysics
D. K. Berry,
1
,
*
M. E. Caplan,
1
,
C. J. Horowitz,
1
,
Greg Huber,
2
,
§
and A. S. Schneider
3
,

1
Center for Exploration of Energy and Matter and Department of Physics, Indiana University, Bloomington, Indiana 47405, USA
2
Kavli Institute for Theoretical Physics, Kohn Hall, University of California, Santa Barbara, California 93106-4030, USA
3
TAPIR, Walter Burke Institute for Theoretical Physics, MC 350-17, California Institute of Technology, Pasadena, California 91125, USA
(Received 28 August 2015; revised manuscript received 6 July 2016; published 1 November 2016)
A striking shape was recently observed for the endoplasmic reticulum, a cellular organelle consisting of
stacked sheets connected by helical ramps [Terasaki
et al.
,
Cell
154
,
285
(
2013
)]. This shape is interesting both
for its biological function, to synthesize proteins using an increased surface area for ribosome factories, and its
geometric 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 magnitude 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 complex 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 conductivities, 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.
DOI:
10.1103/PhysRevC.94.055801
Nuclear pasta, with nucleons arranged into rods, plates, or
other nonspherical shapes, is expected in neutron-star crusts
and in core-collapse supernovae [
1
,
2
]. At just below nuclear
density, these shapes arise from competition between short-
range nuclear attraction and long-range Coulomb repulsion.
Recently we found topological defects consisting of spiral
ramps in molecular dynamics (MD) simulations of nuclear
pasta [
3
,
4
]. Electrons scattering from these spiral ramps could
reduce both the thermal and electrical conductivity of the
neutron-rich crust. This may impact x-ray observations of crust
cooling in transiently accreting neutron stars [
3
], and may also
lead to the decay of neutron-star magnetic fields after about a
million years [
5
].
To explore the energy of these spiral ramp shapes, we start
with the semiempirical mass formula that predicts the binding
energy BE of a nucleus with
A
nucleons and
Z
protons [
6
,
7
],
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 contributions
from the symmetry energy and pairing that will not be
important here. Competition between surface and Coulomb
energy contributions can lead to complex shapes. In this paper,
we study the self-assembly of spiral ramps at a baryon 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
*
dkberry@iu.edu
mecaplan@indiana.edu
horowit@indiana.edu
§
huber@kitp.ucsb.edu

andschn@indiana.edu
fm
3
. Here the system may form flat sheets (lasagne) that
are considerably 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 different
shapes with almost the same energy. For example, if the density
is decreased somewhat, the system may form 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 subdominant terms
in Eq. (
1
). Reinhard
et al.
[
8
], see also [
9
], use a leptodermous
(thin-skinned) expansion and density functionals 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. [
10
]. The curvature energy could
be important for nuclear fission, where a nucleus dramatically
changes shape. However, Reinhard
et al.
find they need to cal-
culate the energy of very large nuclei with
A
in the thousands in
order to extract
a
curve
theoretically. Therefore, this coefficient
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 knowl-
edge) in nuclear physics. We consider a Helfrich-Canham
Hamiltonian
H
0
[
11
,
12
] 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
positive constants representing effective rigidity moduli. This
energy functional was applied in Refs. [
13
,
14
] to biological
2469-9985/2016/94(5)/055801(5)
055801-1
©2016 American Physical Society
BERRY, CAPLAN, HOROWITZ, HUBER, AND SCHNEIDER
PHYSICAL REVIEW C
94
, 055801 (2016)
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, see below, and nuclear pasta. This term may
also stabilize configurations that have additional holes such
as the “nuclear waffle” shapes found in Ref. [
15
]. These
shapes consist of flat sheets with a two-dimensional array
of holes. Alternatively there may be torus or donut-shaped
superheavy nuclei [
16
], 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 spiral
ramp configurations with molecular dynamics simulations
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 and
co-workers in previous works [
3
,
17
21
] and is briefly reviewed
here. It is very similar to a model used by others [
22
,
23
].
Our simulation volume is a cubic box with periodic boundary
conditions which contains pointlike protons 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. Protons and neutrons interact via
the two-body potentials:
V
np
(
r
)
=
ae
r
2
/
+
be
r
2
/
2

,
(3a)
V
nn
(
r
)
=
ae
r
2
/
+
ce
r
2
/
2

,
(3b)
V
pp
(
r
)
=
ae
r
2
/
+
ce
r
2
/
2

+
α
r
e
r/λ
.
(3c)
The
n
and
p
indices indicate whether the potential describes
a neutron-proton, a neutron-neutron, or a proton-proton in-
teraction. Meanwhile,
r
is the separation between each pair
of interacting nucleons,
α
is the fine structure constant, and
quantities
a
=
110 MeV,
b
=−
50 MeV,
c
=−
2 MeV, and

=
1
.
25 fm
2
are parameters of the model. Their values
were chosen in Ref. [
17
] 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 repulsion that causes nuclear
matter to saturate at a density
n
0
=
0
.
16 fm
3
. Finally there is
a long-range (screened) Coulomb repulsion between protons.
Our MD model, Eq. (
3
), predicts a variety of nuclear
pasta shapes such as spheres, rods, or sheets, depending
on for example the density or proton fraction; see Fig. 3
of [
24
]. These shapes are also seen for biological membranes.
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 time step of 2 fm/
c
.
Under these conditions the model tends to form flat sheets, but
these sheets may be connected by ramps.
Similar spiral-shaped ramps appear in cellular biophysics.
Membrane-bound cellular organelles have characteristic
shapes, with the endoplasmic reticulum (ER) being particu-
larly 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 [
13
].
Just as spiral ramps connect the levels of a multilevel “parking
garage,” these so-called “Terasaki ramps” allow sheets to
connect yet remain 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, see Eq. (
2
), largely insensitive to details of
the microscopic interactions [
14
]. In these models of mem-
brane mechanics, the Terasaki ramps are stable, topological
structures, akin to screw dislocations, that affect the entire
morphology of an organelle critical to the metabolism of
eukaryotic cells.
There are, to be sure, dramatic differences between nu-
clear pasta and terrestrial cell biology. Nuclear pasta has a
density near 10
14
g
/
cm
3
, fully 14 orders of magnitude denser
more dense than the aqueous environs of the cell nucleus.
Furthermore, nuclear pasta involves strong interactions be-
tween neutrons and protons in addition to electromagnetic
interactions, while cellular-scale biology is highly screened,
highly overdamped, and dominated by the entropy of water
and complex assemblies of biomolecules. Nonetheless, the
strikingly similar geometry suggests both systems may have
similar coarse-grained dynamics.
To study self-assembly of these ramps in our MD model of
nuclear pasta, we start from initial conditions where the particle
positions are uniformly randomly distributed in the simulation
volume with a Boltzmann velocity distribution. Figure
1
shows the configuration 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 the opacity of the color scale is 0
for
n
p
=
0
.
00 to 0.02 fm
3
and then increases linearly to
1at
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 junctions; 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,
Figs.
1(b)
and
1(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 Figs.
1(e)
and
1(f)
, and the ramps move together
to form the dipole pattern shown in Fig.
2(a)
. This pattern
055801-2
“PARKING-GARAGE” STRUCTURES IN NUCLEAR . . .
PHYSICAL REVIEW C
94
, 055801 (2016)
FIG. 1. Self-assembly of a parking garage structure in a MD
simulation with 40 000 nucleons that started from uniform random
initial positions. Shown in panels (a)–(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.
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
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)
.Thisisalsoa
dipole pattern with the left-handed ramps to the left and the
right-handed ramps to the right. As a result the ramps make
a45
angle with the sheets; see Fig.
2(f)
. Finally a 75 000
particle simulation forms four ramps in a quadrupole pattern
where one left-handed and one right-handed ramp are to the
left as shown in Fig.
2(h)
. For a quadruple pattern the ramps
are observed to make a 90
angle with the sheets; see Fig.
2
(i) and Table
I
.
The arrangement of the helical ramps shown in
Figs.
2(a)
,
2(e)
, and
2(h)
agrees very well with theoretical
predictions in Refs. [
13
,
14
] 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 simulations that are in agreement with
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
FIG. 2. Three axis views,
X
(left),
Y
(center), and
Z
(right), of
the final configuration of three MD simulations. Panels (a)–(c) at top
are for the 40 000 nucleon simulation 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.
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
angle with respect
to the sheets, while ramps in a dipole pattern make an
approximately 45
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 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
.
We compare our results to biological observations. Re-
constructing the three-dimensional (3D) geometry of ER
sheets from two-dimensional serial sections reveals 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 preferred handedness to the ramps, and the observations at
hand, though statistically small, are consistent with right- and
left-handed ramps in equal numbers throughout the organelle.
055801-3
BERRY, CAPLAN, HOROWITZ, HUBER, AND SCHNEIDER
PHYSICAL REVIEW C
94
, 055801 (2016)
FIG. 3. (a) Scanning electron micrograph of a thin slice of
sheetlike ER from a mouse salivary gland. Serial sectioning allows
three-dimensional reconstruction from large numbers of cross-
sections (scale bar
=
200 nm). (b) 3D reconstruction 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 region in (a). (Figure adapted from Terasaki
et al.
[
13
].)
If the sheet edges are treated as effectively one-dimensional
defects (thought to be stabilized 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 minimizing that functional are minimal surfaces,
which have zero mean curvature and locally minimize area.
One consequence of the analysis in Ref. [
14
] 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) minimizes 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 dynamics
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 sheetlike endoplasmic reticulum. Seeing the same helical
shapes in the extraordinarily different systems of nuclear
pasta in neutron stars and in the membranes of eukaryotic
cells strongly suggests that these shapes follow naturally from
common geometric considerations and may be independent of
details of the microscopic interactions.
The curvature energy could help stabilize these helical
ramps, suggesting that they may be common in nuclear pasta.
If that is 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 [
5
]. Alternatively, the helical ramps
that connect flat sheets may allow paired protons to percolate
throughout the system and make it superconducting. In addi-
tion, 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 [
25
]. 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 quasiperiodic
oscillations during magnetar giant flares [
26
]. Finally, the
ramps likely increase the breaking strain, or mechanical
strength, of neutron-star crust. As a result the crust can support
larger mountains that, on a rapidly rotating neutron star,
efficiently radiate gravitational waves [
27
].
One can study soft condensed-matter analogs of nuclear
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 observe 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 biologi-
cally 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 disparate physical
systems allows connections to be made at the deeper level
of symmetry, excitations, and the geometry of topological
defects.
We thank G. Ortiz and S. Setayeshgar for helpful comments.
Part of this work was completed at the Aspen Center for
Physics, which is supported by National Science Foundation
Grant No. PHY-1066293. This research was supported in
part by DOE Grants No. DE-FG02-87ER40365 (Indiana
University) and No. DE-SC0008808 (NUCLEI SciDAC Col-
laboration), and NSF Grant No. PHY11-25915 (KITP, UCSB).
Computer time was provided by the INCITE program. This re-
search used resources of the Oak Ridge Leadership Computing
Facility located at Oak Ridge National Laboratory, which was
supported by the Office of Science of the Department of Energy
under Contract No. DEAC05-00OR22725.
[1] D. G. Ravenhall, C. J. Pethick, and J. R. Wilson,
Phys. Rev. Lett.
50
,
2066
(
1983
).
[2] M. Hashimoto, H. Seki, and M. Yamada,
Prog. Theor. Phys.
71
,
320
(
1984
).
[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] A. S. Schneider, D. K. Berry, M. E. Caplan, C. J. Horowitz, and
Z. Lin,
Phys. Rev. C
93
,
065806
(
2016
).
[5]J.A.Pons,D.Vigano,andN.Rea,
Nat. Phys.
9
,
431
(
2013
).
[6]C.F.vonWeizs
̈
acker,
Z. Phys.
96
,
431
(
1935
).
[7] P. Moller, W. D. Myers, W. J. Swiatecki, and J. Treiner,
At. Data
Nucl. Data Tables
39
,
225
(
1988
).
055801-4
“PARKING-GARAGE” STRUCTURES IN NUCLEAR . . .
PHYSICAL REVIEW C
94
, 055801 (2016)
[8] P.-G. Reinhard, M. Bender, W. Nazarewicz, and T. Vertse,
Phys. Rev. C
73
,
014309
(
2006
).
[9] M. Durand, P. Schuck, and X. Vifias,
Z. Phys. A
346
,
87
(
1993
).
[10] K. Nakazato, K. Iida, and K. Oyamatsu,
Phys. Rev. C
83
,
065811
(
2011
).
[11] W. Helfrich, Z. Naturforsch.
28c
, 693 (1973).
[12] R. Canham,
J. Theor. Biol.
26
,
61
(
1970
).
[13] M. Terasaki
et al.
,
Cell
154
,
285
(
2013
).
[14] J. Guven, G. Huber, and D. M. Valencia,
Phys. Rev. Lett.
113
,
188101
(
2014
).
[15]A.S.Schneider,D.K.Berry,C.M.Briggs,M.E.Caplan,and
C. J. Horowitz,
Phys.Rev.C
90
,
055805
(
2014
).
[16] W. Nazarewicz, M. Bender, S. Cwiok, P. H. Heenen, A. T.
Kruppa, P.-G. Reinhard, and T. Vertse,
Nucl. Phys. A
701
,
165
(
2002
).
[17] C. J. Horowitz, M. A. Perez-Garcia, and J. Piekarewicz,
Phys. Rev. C
69
,
045804
(
2004
).
[18] C. J. Horowitz, M. A. Perez-Garcia, J. Carriere, D. K. Berry,
and J. Piekarewicz,
Phys. Rev. C
70
,
065806
(
2004
).
[19] C. J. Horowitz, M. A. Perez-Garcia, D. K. Berry, and J.
Piekarewicz,
Phys.Rev.C
72
,
035801
(
2005
).
[20] C. J. Horowitz and D. K. Berry,
Phys.Rev.C
78
,
035806
(
2008
).
[21] A. S. Schneider, C. J. Horowitz, J. Hughto, and D. K. Berry,
Phys.Rev.C
88
,
065807
(
2013
).
[22] C. O. Dorso, P. A. Gimenez Molinelli, and J. A. Lopez,
Phys.Rev.C
86
,
055805
(
2012
).
[23] P. N. Alcain, P. A. Gimenez Molinelli, J. I. Nichols, and C. O.
Dorso,
Phys.Rev.C
89
,
055801
(
2014
).
[24] M. E. Caplan and C. J. Horowitz,
arXiv:1606.03646
.
[25] N. Andersson and K. D. Kokkotas,
Int. J. Mod. Phys. D
10
,
381
(
2001
).
[26] A. L. Watts and T. E. Strohmayer,
Adv. Space Res.
40
,
1446
(
2007
).
[27] C. J. Horowitz and Kai Kadau,
Phys. Rev. Lett.
102
,
191102
(
2009
).
055801-5