# Learning to self-fold at a bifurcation

Chukwunonso Arinze1, Menachem Stern2, Sidney R. Nagel1, Arvind Murugan1
Department of Physics, University of Chicago, Chicago, IL 60637  Department of
Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104

###### Abstract

Disordered mechanical systems can deform along a network of pathways that
branch and recombine at special configurations called bifurcation points.
Multiple pathways are accessible from these bifurcation points; consequently,
computer-aided design algorithms have been sought to achieve a specific
structure of pathways at bifurcations by rationally designing the geometry and
material properties of these systems. Here, we explore an alternative physical
training framework in which the topology of folding pathways in a disordered
sheet is changed in a desired manner due to changes in crease stiffnesses
induced by prior folding. We study the quality and robustness of such training
for different ‘learning rules’, that is, different quantitative ways in which
local strain changes the local folding stiffness. We experimentally
demonstrate these ideas using sheets with epoxy-filled creases whose
stiffnesses change due to folding before the epoxy sets. Our work shows how
specific forms of plasticity in materials enable them to learn non-linear
behaviors through their prior deformation history in a robust manner.

Metamaterials [1, 2] and smart materials [3] are often designed to show
specific behaviors. For example, mechanical topological insulators localize
response to forces [4, 5] while elastic networks with allostery communicate
deformations over a long range [6, 7, 8]. More complex mechanical structures
can show a specific deformation, e.g., a smile-shaped deformation, in response
to a spatially textured pattern of forces [9]. Most commonly, mechanical
systems are rationally designed to show such behaviors by searching parameter
space on a computer [10, 11]. However, an alternate approach explored recently
is that of _physical learning_ [12, 13, 14, 15, 16]: during a period of
training, the material is shown examples of the desired behavior, prompting
autonomous changes in the material parameters that promote the desired
behavior. No computers are involved in optimizing the properties of such a
system [17, 18, 19].

Physical learning is a more constrained way of exploring parameter space than
optimization on a computer. However, an autonomous physical learning process
offers the advantage of learning from real examples of stimuli and response
(as opposed to a theoretical specification of the problem) and does not rely
on an accurate model of the material [20]. Physical learning might also allow
for continual learning of new functionalities _in situ_ , as requirements
change [14]. A major open question is what kinds of local adaptive processes
available in real systems constrain physical learning [12]. Recent work has
shown that natural processes within an EVA foam network can train the network
for an auxetic response by simply aging the material in different
configurations [13]. Broader questions remain - what is the impact of
different kinds of local learning dynamics on the feasibility and quality of
learning?

Here, we explore how the quality of physical learning depends on local
adaptive processes (which we call ‘learning rules’) through theoretical
analyses and an experimental demonstration. We focus on training a
fundamentally non-linear feature in marginal mechanical systems, a bifurcation
point [21, 22, 23, 24, 25, 26]. Mechanical bifurcations occur at degenerate
configurations from which multiple non-linear zero modes are accessible.
Bifurcations cannot be described by a linear approximation even for small
deformations since the energy vanishes to two leading orders; they correspond
to singularities of the energy function with an excess of linearized zero
modes and self-stress states that disappear at next to leading order.

Bifurcation points can potentially be exploited to create multi-functional
systems [15] and have been studied in the context of mechanical linkages [27,
28, 29] for robotics and topological meta-materials [25]. However,
bifurcations can also make the system hard to control [30]. For example,
folding outcomes at these singularities can be unpredictable and depend on the
precise spatial pattern of forces used [26, 31, 32, 33]. While generic
mechanical systems show bifurcations in some parts of configuration space [34,
35, 36], bifurcations are especially a problem for thin creased sheets because
the flat state configuration is necessarily degenerate [26, 31]. In
particular, creased thin sheets with nominally 1 degree of freedom (according
to Maxwell counting), often called ‘self-folding origami’, always have a
bifurcation at the flat state that is the meeting point of an exponentially
large number of distinct folding modes. Consequently, such ‘self-folding’
sheets are hard to control at the flat state (despite the name) since the
precise spatial pattern of forces applied will determine the folding geometry.

We focus on the training of such creased thin sheets [37, 38, 39] to
manipulate the bifurcation point in configuration space (Fig.1A); in
particular, we seek to increase the fraction of all force patterns that result
in one folding pathway (i.e., attractor size, Fig.1B) through a physical
learning rule - crease softening due to folding (Fig.1C, Fig.1D). We find that
successful training relies on creating heterogeneity in stiffnesses across the
sheet. However, this heterogeneity can rapidly diminish with further training
for some classes of learning rules while other classes of learning rules,
threshold-like in strain, are more robust. We test some of these ideas with an
experimental prototype in which a sheet with epoxy-filled creases is folded
back and forth along one pathway at a bifurcation. Such folding during the
‘training period’ (i.e., before the epoxy sets) extrudes epoxy to different
extents in different creases, creating a heterogeneous system. We test the
trained sheet by applying different forces and find successful training in
systems with 4 and 7 creases.

## Results

### Training and Stiffness heterogeneity

We begin with a theoretical study of crease patterns made of 4-valent vertices
as shown in Fig.1A. Maxwell constraint counting gives this system 1 degree of
freedom but in reality, this structure has two non-linear 1 degree of freedom
pathways that meet at a bifurcation at the flat state [40]. As a toy model of
energy near a generic bifurcation, consider $E(x,y)=\lambda x^{2}y^{2}$ where
$x,y$ parameterize deformations of a mechanical system. Motions along $x=0$
and $y=0$ are true zero energy pathways, staying at zero energy for large
deformations. However, a linearized analysis at $(x,y)=(0,0)$ suggests a
2-dimensional vector space of zero modes (along with a self-stress state if we
had a mechanistic model); further, any small deformation along $x=0$ or $y=0$
will reduce the zero mode space to 1 dimension (and eliminate the self-stress
mode). Thus, any linearized analysis will fail to identify the true zero modes
$x=0$ and $y=0$. On the other hand, any spatial pattern of forces applied to a
4-vertex [26] at the flat state bifurcation will result in folding along one
of the two zero energy folding pathways.

The non-linear force-response relationship of disordered mechanical systems at
such bifurcations can be summarized by an ‘attractor diagram’, shown
schematically in Fig.1B. The space is a 2-d schematic representation of high-
dimensional space of spatial force patterns. Each colored region represents
one of the discrete folding outcomes that is realized for all force patterns
represented by that region. Earlier work has shown that this attractor
structure can be changed by changing the geometry [26, 31], pre-biasing
preferred directions [41, 42, 43] of the sheet, or controlling the relative
stiffness of different creases [44].

Typically, solving the inverse problem for attractors requires a complex
computer algorithm; for example, Linear or Quadratic Programming algorithms on
a computer [44] can determine crease stiffnesses that eliminate a chosen
pathway in a saddle-node bifurcation. Here, we investigate whether the inverse
problem can be solved by the sheet itself through a natural physical process,
with no computers involved.

Figure 1: The topology of folding pathways at a bifurcation can be changed
through a physical training protocol. (A) Thin creased sheets can fold along
distinct folding pathways that bifurcate from the flat state with outcomes
determined by the spatial force pattern. (B) A 2-d schematic representation of
high-dimensional attractors in the space of force patterns. Each colored
region represents one of the discrete folding outcomes that is realized for
all force patterns represented by that region. Specific attractors can be
enlarged, shrunk or eliminated by changing crease stiffnesses through design
or physical learning. (C,D) We model a physical learning process that changes
a crease $i$’s stiffness $\kappa_{i}$ (about the flat state) based on its
folding strain $\rho_{i}$.

Throughout this paper, each crease, $i$, has a ‘crease stiffness’,
$\kappa_{i}$ which refers to the stiffness for the folding strain (or folding
angle) about the flat state (i.e., unstrained/unfolded state) with rest angle
maintained at zero; i.e., each crease, $i$ has a folding energy
$(1/2)\kappa_{i}\rho_{i}^{2}$ for folding strain $\rho_{i}$ (see Fig.1D) and
crease stiffness $\kappa_{i}$. Thus, this stiffness does not prefer any
folding orientation ($\rho\to-\rho$, which is known as mountain-valley
symmetry) and only depends on the magnitude of the folding strain/angle of the
creases, $\rho_{i}^{2}$. In contrast, a creased sheet of paper develops non-
zero rest angles at each crease that breaks the mountain-valley symmetry
($\rho\to-\rho$) at creases and maintains the sheet away from the flat state
bifurcation, without actually reshaping the bifurcation itself [44]. In this
work, we focus on the non-trivial problem of shaping bifurcations, while
maintaining mountain-valley symmetry at each crease, so the trained sheet can
still be completely flattened. Consequently, our trained sheet can
symmetrically access either the positive or negative components of a given
folding pathway.

Figure 2: Eliminating a select folding pathway for a 4-vertex through physical
training. (A) A single 4-vertex has two distinct folding modes (blue, pink).
(B) The blue pathway’s attractor size during a training process in which the
4-vertex is folded by a fixed force pattern repeatedly. The stiffness
$\kappa_{i}$ of each crease, $i$, are represented by four horizontal lines at
each time point during training; all $\kappa_{i}$ start at the same stiffness
and decrease according to Eq.1. $\kappa_{i}$ are in units of the bulk modulus
of the stiff faces. Crease stiffnesses, $\kappa_{i}$, become heterogeneous
during training, before becoming homogeneous again upon over-training. (C)
Energy landscape of a sheet at different points during training ($\phi$ is an
angular coordinate in 2-dimensional null space; see SI for details). Landscape
before training (circles) and after over-training (pentagons) show two
distinct minima; but when attractor size $\sim$ 1 (squares), pink minimum is
eliminated via a saddle-node bifurcation. Note in the under-trained regime
(triangle), the attractor size is below 1, but well above the initial value of
0.56

As a first pass, we considered learning rules of the type shown in Fig.1C and
Fig.1D that softens different creases based on the current folding strain:

$d\kappa_{i}/dt=-\lambda\rho_{i}^{2}\kappa_{i}.$ (1)

where $\lambda$, the ‘learning rate’, sets the learning timescale.

As a case study, we begin by applying the above learning rule to 4-vertex
shown in Fig.2A; the 4-vertex generically has two distinct folding branches
that meet at the flat state bifurcation[40]. We simulated a training process
in which the vertex initially has creases of equal stiffness $\kappa_{i}^{0}$;
the sheet is folded repeatedly (using a spatial pattern of forces) along the
positive and negative components of one of the two pathways at the
bifurcation. The sheet is folded to a configuration of finite strain and the
parameters $\kappa_{i}$ are updated according to learning rule above for a
time interval $\tau$. The sheet is then relaxed back to the flat state and re-
folded with the negative of the same pattern of forces and the learning
process is continued. See SI for parameters. We test the attractor size of
each pathway throughout training; the attractor size is determined by applying
a library of $50$ randomly chosen spatial pattern of folding torques to the
creases and counting the fraction of folding torques that result in the chosen
folding pathway. We assume that testing does not cause further changes in the
stiffnesses $\kappa_{i}$ in our simulations (though a real system would
continue changing due to testing).

We find that during training, different creases fold to different extents
$\rho_{i}$ in specific ratios characteristic of the branch we fold along.
Consequently, the learning rule Eq.1 creates heterogeneity in the initially
homogeneous creases stiffness $\kappa_{i}$. When crease stiffnesses are
relatively heterogeneous, the attractor size of the chosen pathway increases
from an initial value of $0.56$ to nearly $1$, i.e., nearly all spatial
patterns of forces result in the chosen folding pathway; see Fig.2B. However,
further training reduces the heterogeneity in creases as all stiffnesses
approach zero. In this ‘overtrained’ regime, the attractor size of the desired
pathway drops down to $0.56$ again, no bigger than for the initial untrained
sheet; see Fig.2B.

To illustrate the phenomenon of saddle-node bifurcation, we computed the
energy of the origami sheet in different folded configurations, while its
crease stiffnesses,$\kappa_{i}$ were evolved by the learning rule, Eq.1. We
studied folded configurations within the null space of the origami’s potential
at the flat state. This null space is a 2-dimensional space [44]. We
parameterize this null space with the variables $r$ and $\phi$. We selected a
circle (defined as $r=0.5$) in this null space and computed the energies of
various folded configurations represented by points on the circumference of
this circle. An energy minimum in this energy plot (Fig.2C) corresponds to the
existence of a stable folded configuration. For untrained sheets, we see two
minima corresponding to the two true non-linear folding modes. During
training, the attractor size of one of these true non-linear folding modes is
reduced to zero, as the mode is destroyed in a saddle-node bifurcation. This
is illustrated in Fig.2C (see pink line), as the local minimum on the left
gradually transforms to a local maximum in the course of training. In the
overtrained regime, the eliminated mode (local minimum on the left) reappears
and we find that both pathways are accessible in response to some spatial
patterns of forces (See images of the two modes atop Fig.2C).

Hence, we find that training can solve an inverse problem for non-linear
behavior, namely that of eliminating one select branch in a saddle-node
bifurcation and, thus, changing the topology of pathways. Further, successful
training is correlated with the development of stiffness heterogeneity (see
Fig.2B); this observation is an example of a larger principle that disordered
systems can learn because the information must be stored in the trained
degrees of freedom; homogeneous creases cannot store such information.
Finally, we find that this particular training rule is prone to over-training
and homogenization of creases if training is carried on for too long (see
Fig.2B).

We have found that heterogeneity in the crease stiffnesses stores the learned
information about the desired pathway; but the learning rule that created the
desired heterogeneity also erases that information upon further training.
Similar erasure of trained response was observed in systems like cyclically
sheared Brownian suspensions and charge-density wave conductors [45].

Figure 3: Quality and robustness of learning depend on the learning rule.
(A,B) For the sheet with 13 creases shown, we model a family of ‘training
rules’ with crease stiffness $\kappa_{i}$ that change with strain $\rho_{i}$
(which has been normalized and, thus,ranges from 0 to 1) in a polynomial
$f(\rho_{i})\sim\rho_{i}^{N}$ or threshold
$f(\rho_{i})\sim\rho_{i}^{M}/(\rho_{i}^{M}+T^{M})$ manner. (C) Learning
quality (peak attractor size; see D,F) and robustness (length of training time
over which attractor size $>0.80$; see D,F) for different learning rules are
explored here. (D) With $f(\rho_{i})\sim\rho_{i}^{N}$ rules, the attractor
size of the desired pathway rises to a peak before falling. The peak attractor
size (quality of learning) increases with $N$. (E) The stiffness of each of
the 13 creases, $i$, of the origami are represented by the 13 horizontal
lines. They all start out with the same stiffness and decrease according to
the different polynomial learning rules indicated by the colors of the line.
$\kappa_{i}$ are in units of the bulk modulus of the stiff faces. The
polynomial rules also show a regime of over-training where stiffness
$\kappa_{i}$ become homogeneous again as all $\kappa_{i}$ go to zero. The
horizontal lines representing the crease stiffnesses for the different
polynomial learning rules are plotted at select times, so as to avoid the
lines of the different polynomial learning rules from overlapping. (F)
Attractor size vs training time for threshold learning rules. (G) The
stiffness of each of the 13 creases, $i$, of the origami are represented by
the 13 horizontal lines. They all start out with the same stiffness and
decrease according to the different threshold learning rules indicated by the
colors of the line. $\kappa_{i}$ are in units of the bulk modulus of the stiff
faces. Crease heterogeneity grows with training, and does not decrease (no
overtraining) with continued training for some of the threshold learning rules
(orange line).The horizontal lines representing the crease stiffnesses for the
different threshold learning rules are plotted at select times, so as to avoid
the lines of the different threshold learning rules from overlapping.

### Robustness of different learning rules

Qualitatively, many real materials soften with strain as captured by the
learning rule Eq.1, but might differ quantitatively. Different materials will
have different learning rules; some might soften proportionally to their
strain or to higher powers of their strain and yet others might be even more
sensitive, softening only for strain above a specific threshold. Other non-
equilibrium systems can show more complex learning, where synaptic weights or
learning degrees of freedom can both strengthen or soften over time; we do not
investigate those cases.

We investigated whether such quantitative differences in material properties
might have a qualitative effect on the quality and robustness of learning. We
considered families of rules of the type:

$d\kappa_{i}/dt=-\lambda f(\rho_{i})\kappa_{i}.$ (2)

The first family we considered were different polynomial forms
$f(\rho)\sim\rho^{N}$. Note that Eq.1 is the case where N=2 (See Fig.3A). The
second family, defined by a Hill coefficient $M$,
$f(\rho)\sim\rho^{M}/(\rho^{M}+T^{M})$ corresponds to sigmoidal dependence
often seen in real systems; small strains do not cause significant aging or
change in stiffness but strains above a characteristic threshold $T$ cause
stiffness changes; further the precise amount of strain does not matter beyond
this threshold $T$ (See Fig.3A).

We trained a larger disordered creased sheet with 13 creases and 4 vertices
(See Fig.3B); this sheet had 16 non-linear modes meeting at the flat state
bifurcation. We used the same training protocol as for the single vertex in
Fig.2: we selected one pathway as the desired pathway and applied the learning
rule as the structure was folded repeatedly into the positive and negative
components of the selected pathway. As earlier, we quantified successful
learning by the attractor size of the desired pathway, i.e., folding the sheet
with a library of $50$ random force patterns and computing the fraction of
force patterns that result in a specific folded mode. We assume that
stiffnesses do not change during such testing.

Among polynomial learning rules $f(\rho)\sim\rho^{N}$, we found that training
rules with higher order polynomials (higher $N$) resulted in better training.
We found a higher peak attractor size and training was also more robust -
learning quality stayed higher for longer (See Fig.3C). However, all
polynomial learning rules were still susceptible to overtraining, during which
the crease stiffness heterogeneity was lost (See Fig.3D, Fig.3E). With
threshold-like learning rules, the over-training problem was nearly eliminated
(See Fig.3F). We found that creases that fold less than the threshold $T$ do
not soften at all and hence the learned heterogeneity in stiffness is
maintained over time. However, there is a trade-off; the threshold $T$ of the
learning rule needs to be within the range of strains experienced during
training. If $T$ is too large, no training would occur. If $T$ is too small,
training will fail to create sufficient heterogeneity in stiffnesses
$\kappa_{i}$ to encode information about the desired mode (See Fig.3G).

### Experimental demonstration

With these theoretical results in place, we demonstrated these ideas with an
experimental prototype. While many previous works have implemented creased
sheets in systems ranging from graphene on the nanoscale [46] to mylar and
cardboard on the mesoscale [47] to solar panels on satellites [48], these
works generally have fixed stiffness in different creases and thus an
inability to learn folding behaviors. We note that paper or cardboard is
affected by folding but typically develops a non-zero rest angle upon folding
and is thus likely to fold in the same way again. But as noted earlier, such a
displacement from the bifurcation does not reshape the bifurcation which still
exists if the sheet is forced into the flat state. Here, we create a prototype
that maintains $\rho\to-\rho$ symmetry at each crease and hence can still be
laid flat after training.

We created a sheet with gullies at creases by sandwiching a Tyvek sheet
between acrylic pieces that were laser cut to serve as stiff faces of a crease
pattern. See Fig.4A. Consequently, the creases correspond to gullies of width
$w$ (set by the gap between adjacent acrylic pieces) and depth $h$ (set by the
thickness of acrylic) on both sides of the sheet. A slow setting epoxy
solution is created from a mixture of epoxy resin and a curing agent in the
ratio $1:2$. The creases are filled with epoxy on both sides of the sheet; the
epoxy takes $\sim 60$ minutes to set. See SI for details. During this setting
time (the ‘training phase’), the crease is folded to an angle $+\rho$,
flattened again and folded to angle $-\rho$ in the other direction. Such
folding will extrude an amount of epoxy $h(\rho)$ from the crease gully that
depends on the magnitude of folding $|\rho|$. Only epoxy remaining inside the
crease determines the stiffness of that crease at the end of the training
phase. Thus, the amount of epoxy extruded $h(\rho)$ determines the change in
stiffness $\Delta\kappa$ during training and thus determines the form of the
learning rule in Eq.2. By folding in both directions ($\pm\rho$) during
training, we maintain mountain-valley symmetry and zero rest angles at the
crease; consequently, the trained sheet can still be laid flat.The training
protocol is illustrated in Fig.4B.

We constructed a vertex with 4 creases (studied theoretically in Fig.1),
resulting in two folding modes (pink and blue) shown in Fig.4C. Initially,
before any epoxy is added, all creases are free folding. We tested the
response of this free folding vertex to forces applied at different points
along the boundary of one of the sectors. We found that forces at 3 of 7
locations leads to the blue folding mode while forces at the 4 locations
result in the pink mode (See Fig.4C). Forces applied to all points in the
other sectors led to either the blue or the pink mode. We then filled the
creases with slow setting epoxy, marking the start of the training phase. We
folded the vertex into a selected mode (blue) with characteristic folding
angles $\rho_{i}^{blue}$ at each crease, reverted to the flat state and folded
along the negative component $-\rho_{i}^{blue}$ of the same blue pathway. We
folded until the largest folded crease could not be folded further; in this
way, the magnitude of folding is approximately the same along the positive and
negative components of the pathway and over multiple instances of folding
during training. Throughout this training, the vertex was clamped in a
vertical configuration but was rotated periodically to prevent epoxy from
flowing out of the creases due to gravity. We repeated this folding process
for 60 minutes, folding back and forth along the positive and negative
components of the blue mode. See SI for details of the experiment.

After the epoxy set, we found that forces at all 7 test locations now lead to
the blue mode as illustrated in Fig.4C, showing that the training procedure
had modified the flat state bifurcation, presumably by eliminating the pink
branch at a saddle-node bifurcation away from the flat state [44]. Using a
fresh sample, we also successfully repeated the training process above to
retain the pink mode and eliminate the blue mode instead.

Figure 4: Experimental realization of training through epoxy extrusion. (A)
(schematic) We created creased sheets by sandwiching a thin membrane between
thick acrylic sheets that serve as stiff faces; the resulting creases of width
$w$ and height $h$ are filled with epoxy (orange). Folding by angle $\rho$
before the epoxy sets (the training period) will extrude epoxy; creases with
larger $\rho$ will extrude more epoxy, resulting in lower stiffness $\kappa$
after the epoxy sets. (B) A 4-crease vertex with epoxy filled creases was
trained by repeatedly folding along one of the two pathways during the setting
period; folding is repeated along positive and negative components of the
chosen pathway to avoid any directional folding bias in the creases. (C)
Testing: Folding outcomes are determined for folding forces applied at
different locations. The untrained sheet folds along two distinct pathways
depending on location of applied force (pink and blue dots). The trained sheet
folds along only one pathway for all folding forces (blue dots), demonstrating
an increased attractor for the blue pathway.

To see if the principles behind this simple demonstration are robust enough to
work in more complex disordered systems, we attempted training on a sheet with
7 creases, 2 vertices and thus 4 distinct pathways at the bifurcation; see
Fig.1A. As shown in Fig.5A, the untrained sheet folds along 3 of those 4
pathways for test forces applied to the center of different faces with one
face held clamped. (The fourth pathway requires torques at specific creases
that cannot be realized by forces at a single face in the clamped
configuration we used.) We filled the creases with epoxy and trained with the
same protocol as earlier, folding along a select pathway (the yellow pathway),
flattening the sheet, folding along the negative branch of that pathway and
repeating the process for 60 minutes. After the training process is completed
and the epoxy has set, we tested the sheet with the same test forces applied
to the faces when it was untrained; we now find that all forces lead to
folding along the target yellow pathway as illustrated in Fig.5B. Thus, the
flat state bifurcation has been successfully trained to eliminate the other
pathways, presumably through saddle-node bifurcations away from the flat
state.

Figure 5: Learning to expand a select attractor in a complex sheet. A sheet
with 7 creases and 2 vertices has 4 distinct folding pathways (shown at the
bottom). (A) Before training, three pathways are accessible by forces applied
to different locations shown (circles) or by torques applied to specific
creases (stars). (B) The creases were filled with epoxy and folded back and
forth along the yellow pathway as the epoxy set (the training period). After
training, the sheet was ‘tested’ with the same forces and torques used in (A).
All test forces and torques now result in the yellow pathway, indicating an
expanded attractor size for that pathway.

## Discussion

The study of bifurcations in mechanical systems has attracted attention from
mathematicians [49], roboticists [34, 35] and physicists [31, 26, 44]. Most
work has focused on changing the structure of bifurcations by rationally
changing parameters such as lengths of the elements. Our work here shows at
least some versions of this design problem can be solved by changing stiffness
of joints through a physical training protocol. Our work further suggests that
bifurcations might be physically trainable in mechanical linkages where the
lengths of elements change according to learning rules; changes in length have
been used as a basis for physical training in other contexts [13, 50, 51, 16].

The experimental demonstration here illustrates how a generic physical process
- the extrusion of material at creases - can naturally implement ‘learning
rules’ that confer specific functionality on the system. Other materials
naturally show softening with strain [13], possibly allowing for the
implementation of different functional forms of our learning rules.

The locality of physical learning in mechanical systems contrasts the global
nature of most machine learning algorithm, where learning parameters are non-
locally updated (gradient descent protocols). The training rules are local,
softening each crease in response to strain in that crease. Further, even at
the end of training, the learned stiffness in any one crease does not
immediately favor one folding pathway over another. However, creases that meet
at a vertex [47] have non-linear interactions that constrain their relative
folding; such interacting creases are able to collectively learn and encode
information about a desired pathway even if each individual crease does not
select a pathway by itself.

While we trained for one attractor to grow and occupy most of the force
pattern space, one can also train a system for multiple attractors [15]. Such
training can create a mechanical pattern recognizer, folding into one
configuration in response to one set of force patterns and a distinct
configuration in response to a different set of force patterns. Unlike similar
responsive materials designed on a computer, the learning paradigm here lets
structures learn _in situ_ from real examples of force patterns [13].

###### Acknowledgements.

We thank Heinrich Jaeger, Nidhi Pashine, Chloe Lindeman, Savannah Gowen, and
Martin Falk for discussions. This work was primarily supported by the
University of Chicago Materials Research Science and Engineering Center, which
is funded by National Science Foundation under award number DMR-2011854. M. S.
would like to acknowledge funding from the National Science Foundation via
grant DMR-2005749.

##  I Supplementary Information

###  I.1 Theory

### Theoretical Modeling of Self-Folding Sheets

We model creased sheets using energy-based models developed in previous work
[37, 52, 47]. We assume that creases have a stiffness modeled by torsional
spring elements on each crease [26, 44, 15]. We briefly review the elements of
this model as we build upon this work to simulate the physical learning of
desired folding pathways of self folding origami sheets.

The energy of thin sheet origami is dominated by face bending governed by
mechanical constraints at the origami vertices. Each vertex contributes $3$
constraints on the folding angles of creases around it [52]. Take a vertex
surrounded by $N$ creases denoted with an index $i$ and each folded to an
angle $\rho_{i}$. At the flat state, all $\rho_{i}=0$, which trivially
satisfies all mechanical constraints. One can write down an expansion for
these $3$ non-linear constraints $T_{a}(\rho_{i})$ [26]. The energy of the
folded origami is taken as the sum of squares of the residues of these
constraints $E_{\text{Vertex}}\sim\sum_{a}T_{a}^{2}$, which is independently
summed over different vertices [44]. The energy due to the stiffness in the
creases is $E_{\text{Crease},i}=\frac{1}{2}\kappa_{i}\rho_{i}^{2}$ as
discussed in the main text. Thus, the total energy of a folded sheet is the
sum of vertex and stiff crease energies

$\displaystyle E_{\text{sheet}}(\rho_{i})=\sum_{v\in
vertices}\sum_{a=1}^{3}T_{va}^{2}+\frac{1}{2}\sum_{i\in\text{creases}}\kappa_{i}\rho_{i}^{2}.$
(3)

In the learning protocol presented in this work, the crease geometry is fixed
and so are the vertex constraints $T_{va}$. The change in the energy of a
folded configuration $\rho_{i}$ during training stems directly from the change
in the individual crease stiffness values

$\displaystyle\frac{dE}{dt}=\frac{\partial
E}{\partial\kappa_{i}}\frac{d\kappa_{i}}{dt}=\frac{1}{2}\rho_{i}^{2}\frac{d\kappa_{i}}{dt}.$
(4)

### Simulated Sheet Folding

Using the energy model described previously, we simulated the folding of the
self folding origami via several numerical folding methods [44].

Origami sheets are numerically folded similarly to the way described in [15].
The folded sheet’s configurations are initialized at

$\rho_{i0}=\rho\frac{\tau_{i}}{||\tau_{i}||}$

by a set of external folding torques ($\tau_{i}$) on the creases with
$\rho\equiv||\rho_{i}||$ the folding magnitude, chosen to be $\rho=0.5$.
However, this initialization point is typically not an energy minimum on the
surface of the hyper-sphere of radius $\rho$. Thus, we numerically relax the
sheet to a local minimum of Eq.3 using MATLAB, subject to a constraint that
fixes the folding magnitude $\rho$:

$\displaystyle\begin{aligned}
&\underset{\rho_{i}}{\text{minimize}}&&E_{\text{sheet}}(\rho_{i})\\\
&\text{subject to}&&||\rho_{i}||=\rho.\end{aligned}$ (5)

This protocol mimics the experimental fast folding of origami sheets, and the
clamping of one crease at a specific folded dihedral angle. It was tested and
validated in [44]. The results of this folding protocol are similar to torque
based folding of sheets using Newtonian methods.

### Simulation of Sculpting Folding Pathways Through Physical Learning

Starting with Eq.4 above we see that the evolution of the energy landscape
and, thus, folding pathway is driven by $\frac{d\kappa_{i}}{dt}$. Physical
learning is introduced by the dynamic specification of
$\frac{d\kappa_{i}}{dt}$ defined by various physical learning rules, which are
functions of the fold angles of the mode $\rho^{Teacher}_{i}$, whose attractor
size we desire to expand:

$\displaystyle d\kappa_{i}/dt=-\lambda f(|\rho^{Teacher}_{i}|)\kappa_{i}.$ (6)

where $\lambda$, the ‘learning rate’, sets the learning timescale and
$\rho^{Teacher}_{i}$ is a vector defining the fold angles of the creases of
the desired mode whose attractor size we want to increase. Note that
$\rho^{Teacher}_{i}$ is obtained by folding the creases of the origami with an
external torque, $F^{Teacher}_{i}$ as described above. We generated the
components of $F^{Teacher}_{i}$ by first randomly selecting a number from a
normal distribution. Next we normalized this vector and used it to fold the
creases of the origami as described above. We checked if the scalar product
between the normalized vector and the resulting normalized folded mode is
greater than 0.99. If it is not, we generate another $F^{Teacher}_{i}$ by
selecting another random set of numbers for its components and check if the
new vector and its resulting folded mode have a scalar product greater than
0.99. If the scalar product is greater than 0.99, then the final resulting
folded mode is then normalized and assigned to $F^{Teacher}_{i}$. The learning
rule is specified by $f(\rho^{Teacher}_{i})$ which can be a linear, quadratic
or a threshold function of $\rho^{Teacher}_{i}$. We simulated a physical
learning process in which the creases initially had a uniform stiffness
$\kappa_{i}^{0}=0.02$ (a unit of stiffness represents the bending modulus of
the stiff faces of the origami) and evolved with a learning rate
$\lambda=0.01$ per training round. Thus, we have specified not just a teacher
for the physical learning process, $F^{Teacher}_{i}$, but also a curriculum or
learning rule, $f(\rho^{Teacher}_{i})$, for the physical learning process.

### Simulation of Testing Protocol

After each round of training via physical learning as described above, the
attractor size of the desired mode is computed. To calculate the attractor
size, a set of an array of test torques, $TestTorqueSet$ is created. Each
element of the set, $TestTorqueSet$ is defined by a vector $F^{external}_{i}$,
whose components, $F^{external}_{i,j}$, represent the magnitude of the torque
applied to each crease. Each vector $F^{external}_{i}$ is normalized and used
to fold the creases of the self folding origami sheet as described above.
Folding with each external torque, $F^{External}_{i}$, results in a folded
mode $\rho^{Folded}_{i}$. Note that the folded mode vector,
$\rho^{Folded}_{i}$ is also normalized as well. This folded mode,
$\rho^{Folded}_{i}$ is compared to the desired mode, $\rho^{Teacher}_{i}$
whose attractor size we seek to increase.

To compare the folded mode, $\rho^{Folded}_{i}$, with the desired mode,
$\rho^{Teacher}_{i}$, we take the scalar product between both vectors. If the
scalar product between the two vectors is above 0.9, then we consider the two
modes as similar, and one and the same. We count the number of external
torques, $F^{External}_{i}$, in the set $TestTorqueSet$, whose folded mode,
$\rho^{Folded}_{i}$, are considered similar to $\rho^{Teacher}_{i}$. We then
express this count number as a fraction of the cardinality of the set
$TestTorqueSet$. This fraction defines the attractor size of the desired mode.

### Quality of Learning

Two agents drive the physical learning process: the teacher,
$F^{Teacher}_{i}$, and the curriculum/learning rule, $f(\rho^{Teacher}_{i})$.
We found the kind of teacher selected does not affect the quality of the
physical learning as long as it results in the desired folded mode. Hence, the
quality of learning is determined by the kind of learning rule selected. We
quantify the quality of learning for various learning rules with two
parameters: the peak attractor size attainable and the time robustness of the
learning rule. The peak attractor size compares the maximum attractor size
achieved for a desired mode for the different learning rules. Meanwhile, the
time robustness measures the percentage of the training round for which the
self folding origami is optimally trained (i.e. attractor size is above 0.80).
The time robustness is a measure of the training protocol’s resilience against
over-training.

### Calculation of Energy Landscape for A Single Vertex Origami

To illustrate the mechanism by which physical learning alters the energy
landscape of the self folding origami via a saddle-node bifurcation, we
plotted the energy, Eq.3 of the different folded configurations at several
points during the training of the origami. After each round of training, the
stiffnesses in the creases of the single vertex changes and the energy
landscape of the folded configuration space is re-computed.This time-energy
landscape plot shows the elimination of the unwanted folding pathway (mode)
via a saddle-node bifurcation, and the preservation of the desired folding
pathway (mode) after several rounds of training. Further training results in a
recovery of the previously eliminated mode.

###  I.2 Experiments

### Acrylic sheet setup

To create a system naturally capable of learning, we exploited an origami
system with fresh epoxy totally filled into the crease pattern of the origami.
This epoxy is extruded from the creases during the folding of the origami
during the training protocol. This results in a final stiffness (after the
epoxy sets) that depends on the amount of folding of each of the creases. We
laser cut origami patterns in acrylic sheets of thickness 1.5mm; crease lines
were designed to have a gap (or width) of 30mm. Test holes of 10mm diameter
are laser cut on the acrylic sheet at various strategic positions (along the
circumference of the 90∘ plate for the single vertex and on the center of each
plate of the two-vertex). Two copies of such acrylic patterns were each glued
to both sides of a sheet of Tyvek.The corresponding plates of the origami
patterns on each side of the sheet of Tyvek are lined up with each other
before the glue sets. After the glue is set, holes for applying testing forces
on the faces of the origami plates are perforated. The resulting setup has
stiff faces (bending stiffness set by acrylic) and soft creases (stiffness set
by the Tyvek sheet). Origami patterns studied were for a single-vertex and
two-vertex. The single vertex has four creases radiating from a single vertex
at the center of the pattern. The creases of this single vertex pattern form
sector-angles 150∘, 60∘, 90∘, and 60∘. The two-vertex pattern with a total of
7 creases consists of two internal vertices; one vertex has the following
sector angles 107∘, 123∘, 82∘, and 48∘, while the second vertex is surrounded
by sector angles 82∘, 54∘, 99∘, and 125∘. The two vertices are connected by a
common crease. This connecting crease serves as the boundary dividing the 82∘
sector plate from the 48∘ sector plate of the first vertex, and the 82∘ sector
plate from the 54∘ sector plate of the second vertex. The sector angles,
crease lengths, and position of the vertices for both the single-vertex and
two-vertex are specified in a supplementay PDF file which can be used to laser
cut these patterns.

### Epoxy and training

#### Epoxy Mixture

An epoxy solution is made by mixing epoxy resin with its curing agent in the
ratio 1:2. This mixture is stirred for about five minutes and poured into the
creases of the origami pattern on both sides of the assembly. Note, that if
the epoxy had been mixed in the ratio 1:1, upon curing, it would be so stiff
that the origami assembly would be difficult to fold, without destroying the
assembly; such epoxy when hardened is also brittle and would fracture under a
bending moment attempting to fold the origami assembly. Hence, we mix epoxy
resin and curing agent in the ratio 1:2, allowing for crease folding upon
curing of the epoxy mixture, without disintegration of the origami assembly.

#### Folding of Origami Assembly

The origami assembly filled with watery epoxy in its creases is folded into
the desired configuration that is to be trained for. The origami assembly is
trained by folding the assembly back and forth along the positive and negative
components of a desired folding pathway. This cyclic folding between both
pairs of the desired mode is repeated for an hour, during which the epoxy
solution begins to cure and is no longer watery. A simple folding protocol is
utilized to fold the origami assembly into the desired configuration: one of
the plates of the origami assembly is fixed while pushing or pulling on any of
the other plates of the origami assembly with a normal force or a turning
torque exerted at a single contact point on any of the non-fixed origami
plates.

#### Training Under Gravity

Since the epoxy is still watery during origami training, it needs to be
trained on a rotating platform to avoid epoxy from flowing out of the creases
due to the influence of gravity. The rotating platform consists of two
standing laboratory clamps screwed to the optical table and situated 600mm
apart. A rod is horizontally supported by the claws of the two standing
laboratory clamps, but the rod is allowed to freely rotate within the claws of
the clamp. A lab clamp retort is then fixed clamped (allows for no rotation or
slipping) to the rotation-free horizontally placed rod with the claws on one
end of lab clamp retort, while the claws on the other end of the lab clamp
retort is clamped to one of the plates of the origami (the plate fixed during
training as previously described). As one hand is used to exert a normal force
or turning torque on one of the free plates of the origami assembly in order
to fold it into the positive and negative components of the desired mode, the
other hand is used to manually turn the rotation-free horizontal rod. This
combined process results in the folding of the origami assembly, while under
rotation, and, thus, prevents the flow of the epoxy solution from the
origami’s creases during this training; while under the influence of gravity.

### Emergence of Crease Stiffness

After training the origami assembly, the origami samples are allowed to
further cure and are left hanging in the lab for a week, allowing the epoxy
solution in the creases to harden, thus, producing an effective stiffness on
the creases of the origami.

### Testing protocol

The plate of the origami assembly fixed during training is clamped. A 100mm
long thread knotted on one end is passed through each of the testing holes of
the origami assembly. The threads are gently pulled normal to the surface of
the origami plate. Upon pulling each thread in each hole, the origami folds
into a given configuration. The resulting configuration for each pull is
recorded. The attractor size of the different folding modes of the origami
assembly is computed. This process is repeated for both the trained and
untrained samples of the origami assembly. The attractor size of the chosen
trained mode before and after training are compared to one another.

## References

  * [1] Yongmin Liu and Xiang Zhang.  Metamaterials: a new frontier of science and technology.  Chemical Society Reviews, 40(5):2494–2507, 2011.
  * [2] Pedro M Reis, Heinrich M Jaeger, and Martin van Hecke.  Designer matter: A perspective.  Extreme Mech. Lett., 5:25–29, 2015.
  * [3] Mukesh V Gandhi and BD Thompson.  Smart materials and structures.  Springer Science & Business Media, 1992.
  * [4] Katia Bertoldi, Vincenzo Vitelli, Johan Christensen, and Martin Van Hecke.  Flexible mechanical metamaterials.  Nature Reviews Materials, 2(11):1–11, 2017.
  * [5] Adrien Saremi and Zeb Rocklin.  Topological elasticity of flexible structures.  Physical Review X, 10(1):011052, 2020.
  * [6] Jason W Rocks, Nidhi Pashine, Irmgard Bischofberger, Carl P Goodrich, Andrea J Liu, and Sidney R Nagel.  Designing allostery-inspired response in mechanical networks.  Proceedings of the National Academy of Sciences, 114(10):2520–2525, 2017.
  * [7] Tsvi Tlusty, Albert Libchaber, and Jean-Pierre Eckmann.  Physical model of the genotype-to-phenotype map of proteins.  Physical Review X, 7(2):021037, 2017.
  * [8] Le Yan, Riccardo Ravasio, Carolina Brito, and Matthieu Wyart.  Architecture and coevolution of allosteric materials.  Proceedings of the National Academy of Sciences, 114(10):2526–2531, 2017.
  * [9] Corentin Coulais, Eial Teomy, Koen De Reus, Yair Shokef, and Martin Van Hecke.  Combinatorial design of textured mechanical metamaterials.  Nature, 535(7613):529–532, 2016.
  * [10] Kenneth Diest.  Numerical methods for metamaterial design, volume 127.  Springer, 2013.
  * [11] Nicholas E Jackson, Michael A Webb, and Juan J de Pablo.  Recent advances in machine learning towards multiscale soft materials design.  Current Opinion in Chemical Engineering, 23:106–114, 2019.
  * [12] Nidhi Pashine.  Local rules for fabricating allosteric networks.  Physical Review Materials, 5(6):065607, 2021.
  * [13] Nidhi Pashine, Daniel Hexner, Andrea J Liu, and Sidney R Nagel.  Directed aging, memory, and nature’s greed.  Science Advances, 5(12):eaax4215, 2019.
  * [14] Menachem Stern, Matthew B Pinson, and Arvind Murugan.  Continual learning of multiple memories in mechanical networks.  Physical Review X, 10(3):031044, 2020.
  * [15] Menachem Stern, Chukwunonso Arinze, Leron Perez, Stephanie E. Palmer, and Arvind Murugan.  Supervised learning through physical changes in a mechanical system.  Proceedings of the National Academy of Sciences, 117(26):14843–14850, 2020.
  * [16] Menachem Stern, Daniel Hexner, Jason W Rocks, and Andrea J Liu.  Supervised learning in physical networks: From machine learning to learning machines.  Physical Review X, 11(2):021045, 2021.
  * [17] Sam Dillavou, Menachem Stern, Andrea J Liu, and Douglas J Durian.  Demonstration of decentralized, physics-driven learning.  arXiv preprint arXiv:2108.00275, 2021.
  * [18] JF Wycoff, S Dillavou, M Stern, AJ Liu, and DJ Durian.  Desynchronous learning in a physics-driven learning network.  The Journal of Chemical Physics, 156(14):144903, 2022.
  * [19] Menachem Stern, Sam Dillavou, Marc Z. Miskin, Douglas J. Durian, and Andrea J. Liu.  Physical learning beyond the quasistatic limit.  Phys. Rev. Research, 4:L022037, May 2022.
  * [20] Logan G Wright, Tatsuhiro Onodera, Martin M Stein, Tianyu Wang, Darren T Schachter, Zoey Hu, and Peter L McMahon.  Deep physical neural networks trained with backpropagation.  Nature, 601(7894):549–555, 2022.
  * [21] Jon Kieffer.  Differential analysis of bifurcations and isolated singularities for robots and mechanisms.  IEEE T. Robotic. Autom., 10(1):1–10, 1994.
  * [22] P Kumar and S Pellegrino.  Computation of kinematic paths and bifurcation points.  International Journal of Solids and Structures, 37(46-47):7003–7027, 2000.
  * [23] Yan Chen and Woon Huei Chai.  Bifurcation of a special line and plane symmetric bricard linkage.  Mech. Mach. Theory, 46(4):515–533, 2011.
  * [24] Davide Bigoni.  Nonlinear solid mechanics: bifurcation theory and material instability.  Cambridge University Press, 2012.
  * [25] Bryan Gin-ge Chen, Nitin Upadhyaya, and Vincenzo Vitelli.  Nonlinear conduction via solitons in a topological mechanical insulator.  Proc. Natl. Acad. Sci. U. S. A., 111(36):13004–13009, 2014.
  * [26] Menachem Stern, Matthew B. Pinson, and Arvind Murugan.  The complexity of folding self-folding origami.  Phys. Rev. X, 7:041070, Dec 2017.
  * [27] David H Myszka, Andrew P Murray, and Charles W Wampler.  Computing the branches, singularity trace, and critical points of single Degree-of-Freedom, Closed-Loop linkages.  J. Mech. Robot., 6(1):011006, February 2014.
  * [28] Lin Li, David H Myszka, Andrew P Murray, and Charles W Wampler.  Using the singularity trace to understand linkage motion characteristics.  In ASME 2013 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, pages V06AT07A063–V06AT07A063. Design Engineering Division, August 2013.
  * [29] D Rocklin, Vincenzo Vitelli, and Xiaoming Mao.  Folding mechanisms at finite temperature.  Preprint at http://arXiv.org/abs/1802.02704, 2018.
  * [30] Harry G Kwatny, Bor-Chin Chang, and Shiu-Ping Wang.  Static bifurcation in mechanical control systems.  In Bifurcation Control, pages 67–81. Springer, 2003.
  * [31] Bryan Gin-ge Chen and Christian D Santangelo.  Branches of triangulated origami near the unfolded state.  Phys. Rev. X, 8(1):011034, 2018.
  * [32] A Gillman, K Fuchi, and PR Buskohl.  Truss-based nonlinear mechanical analysis for origami structures exhibiting bifurcation and limit point instabilities.  International Journal of Solids and Structures, 147:80–93, 2018\.
  * [33] ME Lee-Trimble, Ji-Hwan Kang, Ryan C Hayward, and Christian D Santangelo.  Robust folding of elastic origami.  arXiv preprint arXiv:2109.10989, 2021.
  * [34] Charles W Wampler.  Manipulator inverse kinematic solutions based on vector formulations and damped least-squares methods.  IEEE T. Syst. Man Cyb., 16(1):93–101, 1986.
  * [35] Charles W Wampler, Jonathan D Hauenstein, and Andrew J Sommese.  Mechanism mobility and a local dimension test.  Mechanism and machine theory, 46(9):1193–1206, 2011.
  * [36] David H Myszka, Andrew P Murray, and Charles W Wampler.  Mechanism branches, turning curves, and critical points.  ASME 2012 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, pages 1513–1525, January 2012.
  * [37] Thomas C Hull et al.  Modelling the folding of paper into three dimensions using affine transformations.  Linear Algebra and its applications, 348(1-3):273–282, 2002.
  * [38] Tomohiro Tachi.  Generalization of rigid-foldable quadrilateral-mesh origami.  Journal of the International Association for Shell and Spatial Structures, 50(3):173–179, 2009.
  * [39] Levi H Dudte, Etienne Vouga, Tomohiro Tachi, and Lakshminarayanan Mahadevan.  Programming curvature using origami tessellations.  Nature materials, 15(5):583–588, 2016.
  * [40] Scott Waitukaitis, Rémi Menaut, Bryan Gin-ge Chen, and Martin Van Hecke.  Origami multistability: From single vertices to metasheets.  Physical review letters, 114(5):055503, 2015.
  * [41] Tomohiro Tachi and Thomas C Hull.  Self-foldability of rigid origami.  Journal of Mechanisms and Robotics, 9(2):021008, 2017.
  * [42] Ji-Hwan Kang, Hyunki Kim, Christian D Santangelo, and Ryan C Hayward.  Enabling robust self-folding origami by pre-biasing vertex buckling direction.  Advanced Materials, 31(39):0193006, 2019.
  * [43] Ying Zhou, Carlos M Duque, Christian D Santangelo, and Ryan C Hayward.  Biasing buckling direction in shape-programmable hydrogel sheets with through-thickness gradients.  Advanced Functional Materials, 29(48):1905273, 2019.
  * [44] Menachem Stern, Viraaj Jayaram, and Arvind Murugan.  Shaping the topology of folding pathways in mechanical systems.  Nature communications, 9(1):1–8, 2018.
  * [45] Nathan C Keim, Joseph D Paulsen, Zorana Zeravcic, Srikanth Sastry, and Sidney R Nagel.  Memory formation in matter.  Reviews of Modern Physics, 91(3):035002, 2019.
  * [46] Marc Z Miskin, Kyle J Dorsey, Baris Bircan, Yimo Han, David A Muller, Paul L McEuen, and Itai Cohen.  Graphene-based bimorphs for micron-sized, autonomous origami machines.  Proceedings of the National Academy of Sciences, 115(3):466–470, 2018.
  * [47] Matthew B Pinson, Menachem Stern, Alexandra Carruthers Ferrero, Thomas A Witten, Elizabeth Chen, and Arvind Murugan.  Self-folding origami at any energy scale.  Nat. Commun., 8:15477, 2017.
  * [48] Koryo Miura.  Method of packaging and deployment of large membranes in space.  The Institute of Space and Astronautical Science report, (618):1–9, 1985.
  * [49] Denis Jordan and Marcel Steiner.  Configuration spaces of mechanical linkages.  Discrete & Computational Geometry, 22(2):297–315, 1999.
  * [50] Daniel Hexner, Nidhi Pashine, Andrea J Liu, and Sidney R Nagel.  Effect of directed aging on nonlinear elasticity and memory formation in a material.  Physical Review Research, 2(4):043231, 2020.
  * [51] Daniel Hexner, Andrea J Liu, and Sidney R Nagel.  Periodic training of creeping solids.  Proceedings of the National Academy of Sciences, 117(50):31690–31695, 2020.
  * [52] Tomohiro Tachi.  Geometric considerations for the design of rigid origami structures.  In Proceedings of the International Association for Shell and Spatial Structures (IASS) Symposium, volume 12, pages 458–460. Elsevier Ltd, 2010.