Stability of multi-lamellar lipid tubules in excess water

Tripta Bhatia

In the lyotropic phase of lipids with excess water, multilamellar tubules
(MLTs) grow from defects. A phenomenological model for the stability of MLTs
is developed that is universal and independent of the underlying growth
mechanisms of MLTs. The stability of MLTs implies that they are in hydrostatic
equilibrium and stable as elastic objects that have compression and bending
elasticity. The results show that even with 0.1 atm solvent pressure
differences, the density profile is not significantly altered, thus
determining that the stability is due to the trapped solvent. The results are
of sufficient value in relation to lamellar stability models and may have
implications beyond the described MLT models, especially in other models of
membrane systems.

††footnotetext:  Department of Physical Sciences, Indian Institute of Science
Education and Research Mohali, Sector 81, Knowledge City, Manauli, SAS Nagar,
Punjab, 140306 India.

## Introduction

The curvature of biomembranes and the molecular mechanisms that drive the
generation of curvature are important. 1, 2, 3, 4 Bilayer tubules are observed
in lipid-protein compartments such as giant unilamellar vesicles (GUV), in
vitro. In the GUV system, one of the mechanisms for the generation of
curvature is the compositional asymmetry of the leaflet due to the different
amounts of glycolipid GM1.5 Bilayers can also acquire transbilayer asymmetry
through asymmetric adsorption (or desorption) layers formed by the composition
of the solution, resulting in the generation of membrane shapes in the form of
tubules or buds.6, 7, 3 This article describes MLTs as quasi-equilibrium
structures using elasticity theory. MLTs are made up of only one lipid
1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) and grow from defects. This
implies that the bilayers have the same lipid composition and no asymmetry in
terms of this composition. Fig. 1 shows an epi-fluorescence image of stable
cylindrical MLTs rooted in the lamellar reservoir for which stability
conditions are discussed.

Fig. 1:  Two-dimensional view of multilamellar tubules (MLTs) dispersed in
excess water. The brightest region is the lipid reservoir in which the tubules
are rooted. The image was obtained by epi-fluorescence microscopy. The scale
bar is $20\ \mu$m. Fig. 2:  Growth of tubules. After the addition of water,
the growth of the MLT tip is tracked. The edge of the water front moves at a
speed $\sim 5\ \mu$m/s much faster than the growth speed of MLTs.

## Experimental methods

We have prepared lamellar stacks of DOPC lipid (purchased from Sigma).
Approximately 20-50 $\mu$ l of chloroform lipid solution (containing
approximately 0.2 mole% of membrane dye, 1,2-dihexadecanoyl-sn-
glycero-3-phosphoethanolamine, triethylammonium salt, RhPE purchased from
Molecular Probes) is spread on a glass coverlip. The sample is gently dried in
a nitrogen stream and kept covered inside a desiccator overnight with little
vacuum sufficient to hold the chamber tight for the duration. The coated cover
slip with dried sample was glued to a larger cover slip at the edges using
mica spacers of about 100 $\mu$m thickness. Solvent was introduced between the
coverslips of the sample cells by capillary action. After the dry lamellar
stack is hydrated, water molecules are constantly exchanged between the sample
and the bulk; therefore, the hydration gradient acts as a driving force for
tubule growth, as previously described in detail. 8, 9 Fig. 2 shows tubule
growth after water is added. The edge of the water front moves at a speed
$\sim 5\ \mu$m/s much faster than the growth of the tip of the MLTs rooted in
the reservoir. There are two important results of the previous studies: (i)
the need to obtain a necessary first step for the generation of MLTs is the
presence of defects; and (ii) once such a membrane tubule has formed, the open
chamber needs to be closed to stop all growth and retraction dynamics, thereby
achieving hydrostatic equilibrium.
The open edges were sealed with silicone glue immediately after the solvent
filled the entire gap. We found that in the sealed sample cells, an osmotic
equilibrium is reached between the swollen lamellar stack and the excess water
and the tubes remain stable for a couple of days. Sample cells are observed
under a confocal microscope (LEICA TCS-SP2, He-Ne laser 543 nm) equipped with
a 40x dry objective (0.85 N.A.) having a correction collar. Fig. 3 shows a
two-dimensional confocal cross section of the tip and root of the tubules. The
double arrow indicates the incident laser polarization. In a closed (sealed)
chamber, the diameter of the core of these MLTs is extracted from quantitative
image analysis, namely optimum smoothening.10 Fluorescence confocal polarizing
microscopy (FCPM) observations confirmed that the tubules are multilamellar,
are capped at the end and originate from defects. The bright region in the
image of the tubules is the region where the dye molecules have their
absorption and emission dipole moments oriented parallel to the incident laser
polarization. For a closed (sealed) chamber, water-molecule exchange can occur
only within the sample cell, and an equilibrium is reached for which the
number of water molecules leaving the lamellar stack is exactly balanced by
the number of water molecules entering the stack, and no further growth of
tubules occurs. In this situation, the image reveals that MLTs have a wide
range of $r_{\rm{c}}$ and $r_{\rm{o}}$. Given the earlier observations 9, in
particular the long lifetime of the structures, it is reasonable to consider
these as quasi-equilibrium structures.

## Elasticity of tubules

A phenomenological model is proposed to analyze the stability of simple
cylindrical MLTs with a uniform cross section shown in Fig. 3. The pertinent
experimental observations on which the model is based are listed below.

  1. 1.

After the sample cell is sealed, the tubules remain stable for more than one
day. 9, 5, 8 The model applies to the equilibrium shapes achieved by the MLTs
shown in Fig. 3.

  2. 2.

The tubules are capped at the end (Fig. 3). The membrane fluorescence signal
comes from dye molecules oriented in the bilayer with the absorption and
emission dipole moment always oriented in the plane of the bilayer 9, 8. Thus,
all MLT images exhibit this optical effect, namely photoselection, from which
the membrane morphology of MLTs with a well-defined tip can be defined.

  3. 3.

The tubules originate from defects in the lipid reservoir, shown in Fig. 3. 9

  4. 4.

Tubules have a wide range of $r_{\mathrm{c}}$ and $r_{\mathrm{o}}$. In
general, the examined MLTs have $r_{\rm{o}}$ in the range of 1-30 $\mu$m.
Furthermore, in the stability regime, the structures are investigated with a
z-stack confocal that detects the diameter of the core in the range $88\pm 23$
nm and $6860\pm 50$ nm, respectively, using a special image processing
technique, namely optimal smoothening. 10

Fig. 3:  Two-dimensional confocal cross section of the tip and root of the
tubules. The double arrow indicates the incident laser polarization. FCPM
observations confirmed that the tubules are multilamellar, are capped at the
end, and originate from defects. Top panel: Multilamellar tubule (MLT)
observed by fluorescent confocal polarized microscopy (FCPM) with incident
laser polarization (a) oriented perpendicular to the long axis of the tubule
and (b) parallel to the long axis of the tube. The scale bar is $10\ \mu$m
applied to (a,b). Bottom panel: FCPM image of the root of arbitrary MLTs in
the lipid reservoir, which confirms the multilamellar arrangement. The bright
region of the tubules is the region where the dye molecules have their
absorption and emission dipole moments oriented parallel to the incident laser
polarization.

Given the above observations and, in particular, the long life times of the
structures, it is reasonable to consider these as quasi-equilibrium
structures. Therefore, the theory of elasticity is used to understand the
stability and structure of MLTs.

#### Elasticity of stack of bilayers:

The smectic liquid crystals are one-dimensional “solids” composed of fluid
layers exhibiting quasi-long-range order in the orthogonal direction of the
layers. For small distortions, the elastic free energy (in Cartesian
coordinates) is 11, 12, 13.

$F_{\mathrm{el}}\hskip 1.42262pt=\int\hskip 1.42262pt\Big{[}\frac{B}{2}\hskip
1.42262pt\Big{(}\partial_{z}u\Big{)}^{2}\hskip 1.42262pt+\hskip
1.42262pt\frac{\kappa}{2}\hskip 1.42262ptH^{2}\hskip 1.42262pt+\hskip
1.42262pt\kappa_{\mathrm{G}}\hskip 1.42262ptK\Big{]}\hskip 1.42262ptdx\hskip
1.42262ptdy\hskip 1.42262ptdz\hskip 1.42262pt,$ (1)

where $z$ is the solid direction “like” of the layering, $u\hskip
1.42262pt(x,y,z)$ is the layer displacement field. $(B/2)\hskip
1.42262pt\Big{(}\partial_{z}u\Big{)}^{2}$ is the energy density for
compression (or extension) of the layer spacing with $B$ as the compression
modulus. $(\kappa/2)\hskip 1.42262ptH^{2}$ is the energy density for layer
bending, with the modulus of bending $\kappa$ and the mean curvature
$H=[(1/2)(1/R_{1}\hskip 1.42262pt+\hskip
1.42262pt1/R_{2})=\bigtriangledown_{\perp}^{2}u]$, with the principal radii of
curvature $R_{1}$ and $R_{2}$ and
$\bigtriangledown_{\perp}^{2}=(\partial_{x}^{2}+\partial_{y}^{2})$.
$K=(1/R_{1}R_{2})$ is the Gaussian curvature. The Gaussian curvature term
contributes to the energy only if the system under study has a boundary or
undergoes a change in topology (number of handles). The tubule in Fig. 3 has a
single core, with a uniform core radius $(r_{\rm{c}})$ and an outer radius
$(r_{\rm{o}})$. The edge structure of the tubule is probed by changing the
incident laser polarization by $(\pi/2)$ from Fig. 3a to Fig. 3b, which
confirms that the tubule is capped at the end and that the open edges of the
bilayers are not exposed to the solvent. MLTs grow from defects in the
hydrated lamellar stack 9. The main result of our paper is to develop a
curvature-elastic model for the stability of MLTs, as described in the
following.

#### Elasticity of MLTs:

An MLT is an aggregate of multiple bilayers that are rolled coaxially to form
cylindrical structures characterized by an outer radius $r_{\rm{o}}$ and an
inner core radius $r_{\rm{c}}$ that are dispersed in the bulk water in the
sample cell. The tubules are capped at one end, and the other end is rooted
into the lamellar stack. We have ignored the tubule cap energy, since the
tubules are $(10-100)\ \mu$m long. The elastic energy of the tubular lamella
contributes from the compression of the layer and the curvature of the layer
11, 12 with the direction of the layer changed from $z\rightarrow r$ and the
direction of the in-plane changed from $(x,y)\rightarrow(r,\phi)$, where $r$
is the radial direction, $\phi$ is the polar angle and $z$ lies parallel to
the long axis of the tubules. Equation (1) becomes

$F_{\mathrm{el}}\hskip 1.42262pt=\int\hskip 1.42262pt\Big{[}\frac{B}{2}\hskip
1.42262pt\Big{(}\partial_{r}u\Big{)}^{2}\hskip 1.42262pt+\hskip
1.42262pt\frac{\kappa}{2}\hskip
1.42262pt\left(\partial^{2}_{r}u+\frac{1}{r^{2}}\partial^{2}_{\phi}u+\frac{1}{r}\partial_{r}u\right)^{2}\hskip
1.42262pt\Big{]}\hskip 1.42262ptrdr\hskip 1.42262ptd\phi\hskip
1.42262ptdz\hskip 1.42262pt$ (2)

where the symbols have their usual meaning, $\partial_{r}$ stands for the
differential operator d/d$r$ and
$H=\left(\partial^{2}_{r}u+\frac{1}{r^{2}}\partial^{2}_{\phi}u+\frac{1}{r}\partial_{r}u\right)$
is the mean curvature. The Gaussian curvature term is neglected 14. We retain
the subdominant term in the free elastic energy because the compressive and
curvature stresses must be balanced for stability. We define
$\lambda^{2}=(\kappa/B)$ where $\lambda$ and $B$ have the dimension of length
and pressure, respectively.

If the tube is assumed to have a uniform circular cross section along its
length, then the dependence of $u(r,\phi,z)$, on $\phi$ becomes zero for
hydrostatic equilibrium, i.e; $(\partial_{\phi}u=0)$. Hence, for cylindrically
symmetric MLTs with a uniform cross section, the elastic free energy is;

$\begin{split}{F}_{\rm{el}}&=\int_{r,z,\phi}\
r\frac{B}{2}\left[(\partial_{r}u)^{2}+\lambda^{2}\left(\partial^{2}_{r}u+\frac{1}{r}\partial_{r}u\right)^{2}\right]dr\
d\phi\ dz\\\ &=2\pi\int_{r,z}\ rf[u(r),u^{\prime}(r),u^{\prime\prime}(r)]dr\
dz\end{split}$

where $u^{\prime}(r)=\partial_{r}u,\ u^{\prime\prime}(r)=\partial^{2}_{r}u$
and $F_{\rm{el}}$ have the dimension of energy. We consider all possible
variations in ${F}_{\rm{el}}$ induced by a shift in layer displacement from
$u$ to $(u+\varDelta u)$,

$\begin{split}\delta{F}_{\rm{el}}&=2\pi\int_{r,z}\delta\left(rf[u(r),u^{\prime}(r),u^{\prime\prime}(r)]\right)dr\
dz\\\ &=2\pi\int_{r,z}\left[\frac{\partial(rf)}{\partial u}\delta
u+\frac{\partial(rf)}{\partial u^{\prime}(r)}\delta
u^{\prime}(r)+\cdot\cdot\cdot\right]dr\ dz\end{split}$ (3)

### Condition for stability

The tubules are bounded by water at the top, between the bilayers, and inside
the core. To be in hydrostatic equilibrium, the normal stresses at
$r_{\rm{c}}$ and $r_{\rm{o}}$ must balance the fluid pressure inside and
outside the tubule, respectively, as shown in Fig. 4. If the negative and
positive signs show the direction outward and inward from the bilayer plane,
respectively, then

$\begin{split}\sigma(r_{\rm{i}})=-p_{\rm{i}},\\\ -\
\sigma(r_{\rm{o}})=p_{\rm{o}}\ \end{split}$ (4)

where $\sigma_{\rm{i}}$ and $\sigma_{\rm{o}}$ are the normal stresses in the
innermost and outermost layers that counterbalance the pressures $p_{\rm{i}}$
and $p_{\rm{o}}$ of the solvent in $r_{\rm{i}}$ and $r_{\rm{o}}$,
respectively, as shown in Fig. 4. In equilibrium, the Euler-Lagrange equation
$(\delta F_{\rm{el}}/\delta u)=0$ holds within the bulk of the lamellar region
of the MLT giving,

$\begin{split}\left(\lambda^{2}-r^{2}\right)u^{\prime}(r)&-r\left(\lambda^{2}+r^{2}\right)u^{\prime\prime}(r)\\\
&+r^{2}\lambda^{2}\left(ru^{(4)}(r)+2u^{(3)}(r)\right)=0\end{split}$ (5)

The solution to the Euler-Lagrange equation for $(r>0,\lambda>0)$ is given by

$\begin{split}u(r)&=c_{3}\
\lambda^{2}Y_{0}\left(-\frac{ir}{\lambda}\right)+c_{2}\
\lambda^{2}\left(I_{0}\left(\frac{r}{\lambda}\right)-1\right)+c_{1}\
\log(r)+c_{4}\\\ &=c_{2}\
\lambda^{2}\left(I_{0}\left(\frac{r}{\lambda}\right)-1\right)+c_{1}\log(r)+c_{4},\end{split}$

where $I_{n}(x)$ and $Y_{\rm{n}}(x)$ 15, respectively, represent modified
Bessel functions of the first kind and Bessel functions of the second kind and
of order $n$, and the four integration constants are denoted with the symbol
$c$. The layer displacement $u(r)$ is real. Therefore, we drop the term
proportional to $Y_{0}\left(-ir/\lambda\right)$, which is a pure imaginary.
The constant $c_{4}$ corresponds to the uniform and rigid displacements of the
entire MLT and can therefore be set to zero $(c_{4}=0)$ by an appropriate
choice of the origin of the coordinate system. This gives a certain real value
of $u(r)$.

$u(r)\simeq c_{2}\
\lambda^{2}\left(I_{0}\left(\frac{r}{\lambda}\right)-1\right)+c_{1}\log(r)$
(6)

The gradient of layer displacement or the bilayer compression,
$\partial_{r}u(r)$ is a direct measure of the density change given by

$v(r)=\partial_{r}u(r)=c_{2}\ \lambda\
I_{1}\left(\frac{r}{\lambda}\right)+\frac{c_{1}}{r}$ (7)

Now, it remains to find the two integration constants $c_{1}$ and $c_{2}$
subject to the boundary conditions given by equation (4).

Fig. 4:  An MLT in hydrostatic equilibrium. $\sigma_{\rm{i}}$ and
$\sigma_{\rm{o}}$ are the normal stresses in the innermost and outermost layer
that counterbalance the pressures $\rm{p}_{\rm{i}}$ and $\rm{p}_{\rm{o}}$ of
the solvent. ”W” denotes the water (or solvent), ”L” denotes the lamella,
$r_{\rm{i}}(=r_{\rm{c}})$ and $r_{\rm{o}}$ denote the inner (or core) and
outer radii of the tubule, respectively.

### Solving for $u(r)$

The divergence of normal stress 16 $(\partial_{i}\sigma_{ij})$ on a patch of
membrane with area ($dA=rdrd\phi$) and volume ($dV=rdrd\phi dz$) is related to
the total force per unit volume given by;

$f_{tot}=\int_{V}(\partial_{i}\sigma_{ij})\ dV=\int_{A}\sigma_{ij}\
\hat{n}_{j}\cdot d\hat{A}$ (8)

where $\hat{n}_{j}$ is the normal area vector. Normal stress $\sigma(r)$ and
radial strain $u^{\prime}(r)$ are related to elastic free energy 16 as;

$\begin{split}\sigma(r)&=\frac{\delta{F}_{\rm{el}}}{\delta(u^{\prime}(r))}=\frac{1}{r}\frac{\delta(rf)}{\delta(u^{\prime}(r))}\\\
&=\frac{2\pi
B}{r}\left(\lambda^{2}u^{\prime\prime}(r)+\left(\frac{\lambda^{2}}{r}+r\right)u^{\prime}(r)\right)\\\
&=\frac{2\pi
B}{r}\left(\lambda^{2}v^{\prime}(r)+\left(\frac{\lambda^{2}}{r}+r\right)v(r)\right)\end{split}$
(9)

Solving equations (4) and (9) we obtain the following.

$\begin{split}c_{1}&=r_{\rm{i}}\left[\frac{p_{\rm{i}}r_{\rm{i}}^{2}}{2\pi
B\alpha}+c_{2}\lambda\left(\frac{\lambda
r_{\rm{i}}I_{0}\left(\gamma_{\rm{i}}\right)}{\alpha}-I_{1}\left(\gamma_{\rm{i}}\right)\right)\right],\\\
c_{2}&=\frac{p_{\rm{i}}r_{\rm{i}}^{3}\beta-
p_{\rm{o}}r_{\rm{o}}^{3}\alpha}{2\pi
B\lambda\left[\alpha\left(\beta\left(r_{\rm{i}}I_{1}\left(\gamma_{\rm{i}}\right)-r_{\rm{o}}I_{1}\left(\gamma_{\rm{o}}\right)\right)+\lambda
r_{\rm{o}}^{2}I_{0}\left(\gamma_{\rm{o}}\right)\right)-\lambda
r_{\rm{i}}^{2}\beta I_{0}\left(\gamma_{\rm{i}}\right)\right]}\end{split}$

where $\alpha=\left(2\lambda^{2}+r_{\rm{i}}^{2}\right)$,
$\beta=\left(2\lambda^{2}+r_{\rm{o}}^{2}\right)$,
$\gamma_{\rm{i}}=\left(r_{\rm{i}}/\lambda\right)$, and
$\gamma_{\rm{o}}=\left(r_{\rm{o}}/\lambda\right)$. Substituting $c_{1}$ and
$c_{2}$ into equation (6), $u(r)$ is calculated from equation (6).

### Hydrostatic equilibrium

The bilayer compression is plotted for an MLT with a given pressure difference
$(p_{\rm{i}}-p_{\rm{o}})$ of the solvent. I have chosen the dimensionless unit
of length as $\gamma_{\rm{i}}=(r_{\rm{i}}/\lambda)=1$, and measure the
pressure in units of $B$. For the $L_{\alpha}$ phase under consideration,
$\lambda$ is of the order of a few layer spacings and $B\simeq\ 10$ atm.
The $\partial_{r}u(r)$ for an MLT with the same inner and outer pressures
$(p_{\rm{i}}=p_{\rm{o}})$ of the solvent is plotted and shows that the MLT has
varying non-zero bilayer compression, as shown in Fig. 5 by the black color
curve. If the inner and outer pressures are zero, the bilayer compression,
$\partial_{r}u(r)$ and the stresses are zero throughout the tubule.

Fig. 5: Red curve [$(\rm{p}_{\rm{i}}-\rm{p}_{\rm{o}})=0.1\ $ atm], black curve
($\rm{p}_{\rm{i}}=\rm{p}_{\rm{o}}$). The tubule has $(r_{\rm{c}},\
r_{\rm{o}})=(1,3)$ in dimensionless units of length as $(r_{\rm{i}}/\lambda)$.

In reality, the solvent pressures inside and outside the tubule can be
different, because $(r_{\rm{c}}<r_{\rm{o}})$. A lower $r_{\rm{c}}$ implies a
higher normal stress of the layer curvature (the mean curvature is
1/2$r_{\rm{c}}$ at the core-solvent interface). The normal stress at the outer
radius is clearly lower, as $(r_{\rm{o}}>r_{\rm{c}})$. Thus, it is reasonable
to assume that the solvent pressure inside the tubule (in the core) is greater
than that outside the tubule. We have plotted $\partial_{r}u(r)$ for an MLT
with unequal inner and outer pressures $(p_{\rm{i}}\neq p_{\rm{o}})$ of the
solvent, e.g., $(p_{\rm{i}}-p_{\rm{o}})=0.1\ $ atm. Even in this extreme case,
we find that the compression and, therefore, the density profile of the
lamellar material within the MLT do not change significantly, as shown in Fig.
5, by the red color curve. In such a scenario, the pressure difference can be
stabilized by trapping the solvent within small closed regions that are
dynamically formed in the reservoir.

Fig. 6:  Different size and morphology of multilamellar tubules (MLTs)
dispersed in excess water. The cross-sectional view ($279\ \mu$m on each side)
is obtained by fluorescence confocal polarized microscopy (FCPM) with incident
laser polarization oriented along the vertical axis of the image.

## Conclusions

Stability of MLTs implies that they are in hydrostatic equilibrium and stable
as elastic objects that have compression and bending elasticity. In this
regime, the model is universal and independent of the underlying growth
mechanisms of MLTs. The MLTs did not exhibit any significant motion on the
typical time scales of the experiments in the closed sample cell. Furthermore,
in the stability regime, MLTs are treated as aggregates of lipid bilayers with
an equilibrium spacing of $D_{\rm{eq}}$ between adjacent bilayers, which is
achieved by balancing the intermolecular forces within the stack. In the
presence of defects, growth is driven by the hydration gradient.8 The size and
morphology of MLTs in the experiments cannot be controlled, as shown in Fig.
6. However, if the lamellar stacks are spin coated, then MLTs do not form.9
Multilateral tubule growth, also known as myelin figures (MF) 17, during lipid
hydration is explored using surfactants such as Aerosol-OT, triethylene glycol
monododecyl ether 18, 19, 20, 21, 22, 23, 24, 25, 26, from vesicular
dispersions of lung surfactant extract 27, lipid vesicles with hydrophilic
polymers 28, 29 and aqueous multicomponent polymer solutions 30. In living
systems, an important example is the presence of MF in the pulmonary lining,
an extracellular lipid-protein coat on the alveolar surface, the presence of
which helps reduce surface tension when the lungs are deflated 31, 32, 33. The
model is simple and might be of interest to the field in understanding tubule
stability in a variety of biological processes/applications. The results are
shown in Fig. 5 and show that even with solvent pressure differences of 0.1
atm, the density profile is not significantly altered, thus determining that
the stability is due to the trapped solvent. I also believe that the results
are important enough for the lamellar stability model to be considered for
publication, with implications beyond the MLT models described, especially
when trying to model other membrane systems.

## Acknowledgment

TB acknowledges a Ph.D. thesis research fellowship from the Raman Research
Institute (RRI) , Bangalore, India. TB acknowledges Prof. Yashodhan Hatwalne
(RRI) , Dr. A. Chaudhuri (IISERM), and Prof. J. H. Ipsen (SDU) for useful
discussions.

## References

  * Jarsch _et al._ 2016 I. K. Jarsch, F. Daste and L. Gallop, _J. Cell Biol._ , 2016, 214, 375.
  * Rothman and Lenard 1977 J. E. Rothman and J. Lenard, _Science_ , 1977, 195, 743–753.
  * Lipowsky 2022 R. Lipowsky, _Adv. Biology_ , 2022, 6, 2101020.
  * Bhatia 2022 T. Bhatia, _J. Membr. Biol._ , 2022.
  * Bhatia _et al._ 2018 T. Bhatia, J. A. Canalejo, R. Dimova and R. Lipowsky, _ACS Nano_ , 2018, 12, 4478–4485.
  * Bhatia _et al._ 2020 T. Bhatia, S. Christ, J. Steinkühler, R. Dimova and R. Lipowsky, _Soft Matter_ , 2020, 16, 1246.
  * Lipowsky 2022 R. Lipowsky, _Adv. Colloid Interface Sci._ , 2022, 301, 102613.
  * Bhatia 2021 T. Bhatia, _J. Colloid Interface Sci._ , 2021, 584, 706–713.
  * Bhatia _et al._ 2015 T. Bhatia, Y. Hatwalne and N. V. Madhusudana, _Soft Matter_ , 2015, 11, 5641–5646.
  * Bhatia 2018 T. Bhatia, _Eur. Biophys J._ , 2018, 47, 531–8.
  * de Gennes and Prost 1994 P. G. de Gennes and J. Prost, _The Physics of Liquid Crystals_ , Oxford University Press, Oxford, 1994.
  * Chaikin and Lubensky 1998 P. M. Chaikin and T. C. Lubensky, _Principles of condensed matter physics_ , Oxford University Press, 1998.
  * Helfrich 1973 W. Helfrich, _Z. Naturforsch. and C_ , 1973, 28, 693.
  * Safran 2003 S. Safran, _Statistical thermodynamics of Surfaces and Interfaces and Membranes_ , West View Press, Cambridge, 2003.
  * Abramowitz and Stegun 1972 M. Abramowitz and I. Stegun, _Handbook of Mathematical Functions_ , Dover, New York, 1972.
  * Landau _et al._ 1986 L. D. Landau, L. P. Pitaevskii, E. M. Lifshitz and A. M. Kosevich, _Theory of Elasticity_ , Pergamon Press, 1986.
  * Chapman _et al._ 1966 Chapman, D., Fluck and D. J., _J. Cell Biol._ , 1966, 30, 1.
  * Dave _et al._ 2003 H. Dave, M. Surve and J. Bellare, _J. Colloid Interface Sci_ , 2003, 264, 76–81.
  * Mishima and Yoshiyama 1987 K. Mishima and K. Yoshiyama, _Biochim.Biophys.Acta_ , 1987, 904, 149–153.
  * Huang _et al._ 2005 J. R. Huang, L. N. Zou and T. A. Witten, _Eur. Phys. J. E._ , 2005, 18, 279.
  * Sakurai _et al._ 1989 I. Sakurai, T. Suzuki and S. Sakurai, _Biochimica et Biophysica Acta_ , 1989, 985, 101.
  * Sakurai and Kawamura 1985 I. Sakurai and Y. Kawamura, _Biochimica et Biophysica Acta_ , 1985, 815, 149–152.
  * Zou and Nagel 2006 L. N. Zou and S. R. Nagel, _Phys. Rev. Lett._ , 2006, 96, 138301.
  * Zou 2009 L. N. Zou, _Phys. Rev. E._ , 2009, 79, 061502.
  * Buchanan _et al._ 2000 M. Buchanan, S. U. Egelhaaf and M. E. Cates, _Langmuir_ , 2000, 16, 3718–3726.
  * Reissig _et al._ 2010 L. Reissig, D. J. Fairhurst, J. Leng, M. E. Cates, A. R. Mount and S. U. Egelhaaf, _LANGMUIR_ , 2010, 26, 15192–99.
  * Andersson _et al._ 2018 J. M. . Andersson, K. Roger, M. Larsson and E. Sparr, _ACS Cent. Sci._ , 2018, 4, 1315–1325.
  * Tsafrir _et al._ 2003 I. Tsafrir, Y. Caspi, M. A. Guedeau-Boudeville, T. Arzi and J. Stavans, _Phys. Rev. Lett._ , 2003, 91, 138102.
  * Tsafrir _et al._ 2001 I. Tsafrir, D. Sagi, T. Arzi, M. A. Guedeau-Boudeville, V. Frette, D. Kandel and J. Stavans, _Phys. Rev. Lett._ , 2001, 86, 1138–1141.
  * Preston _et al._ 1980 B. N. Preston, T. C. Laurent, W. D. Comper and G. J. Checkley, _Nature_ , 1980, 287, 499–503.
  * R. J. Sanderson 1977 A. E. V. R. J. Sanderson, _Journal of Cell Biology_ , 1977, 74, 1027–1031.
  * Vangolde _et al._ 1988 L. M. G. Vangolde, J. J. Batenburg and B. Robertson, _Physiological Reviews_ , 1988, 68, 374–455.
  * Weibel _et al._ 1966 E. R. Weibel, G. S. Kistler and G. Tondury, _Zeitschrift Fur Zellforschung Und Mikroskopische Anatomie_ , 1966, 69, 418–427.


# Wavelet-based Loss for High-frequency Interface Dynamics

Lukas Prantl
Technical University of Munich
<EMAIL_ADDRESS>
& Jan Bender
RWTH Aachen University
Tassilo Kugelstadt
RWTH Aachen University
& Nils Thuerey
Technical University of Munich

###### Abstract

Generating highly detailed, complex data is a long-standing and frequently
considered problem in the machine learning field. However, developing detail-
aware generators remains an challenging and open problem. Generative
adversarial networks are the basis of many state-of-the-art methods. However,
they introduce a second network to be trained as a loss function, making the
interpretation of the learned functions much more difficult. As an
alternative, we present a new method based on a wavelet loss formulation,
which remains transparent in terms of what is optimized. The wavelet-based
loss function is used to overcome the limitations of conventional distance
metrics, such as L1 or L2 distances, when it comes to generate data with high-
frequency details. We show that our method can successfully reconstruct high-
frequency details in an illustrative synthetic test case. Additionally, we
evaluate the performance when applied to more complex surfaces based on
physical simulations. Taking a roughly approximated simulation as input, our
method infers corresponding spatial details while taking into account how they
evolve. We consider this problem in terms of spatial and temporal frequencies,
and leverage generative networks trained with our wavelet loss to learn the
desired spatio-temporal signal for the surface dynamics. We test the
capabilities of our method with a set of synthetic wave function tests and
complex 2D and 3D dynamics of elasto-plastic materials.

##  1 Introduction

In the context of generative neural networks, simple distance metrics such as
L1 or L2 distances play an important role. They often serve as a loss function
to validate the output in relation to the ground truth. However, these
conventional distance metrics usually act like a low-pass filter. Small high-
frequency details are thus weighted too weakly and the results become smooth
and can sometimes be very different from sharp, detailed ground-truth data.
This is especially true for multi-modal problems that are ill-conditioned and
where features cannot be unambiguously reconstructed from the information-poor
input. Representative for this are super-resolution tasks, where high-
resolution data are to be reconstructed from low-resolution data (Johnson et
al., 2016; Dosovitskiy & Brox, 2016; Bruna et al., 2016). To counteract the
limitations of conventional distance functions, we introduce a novel wavelet-
based loss function.

Our central motivation is to put more emphasis on high-frequency features. For
the proposed wavelet loss, we therefore first transform the data into wavelet
space and then scale it such that the focus lies on the high-frequency
content. This successfully negates the low-pass effect of the simpler distance
metrics. We choose the wavelet transform (Daubechies, 1992) because it better
fits the nature of many natural data sets, where high-frequency details do not
appear uniformly, but instead are sporadic and isolated. The first choice for
frequency analysis would usually be a Fourier transform, but this usually
expects periodic signals whose frequency does not change significantly over
time. An adaptation of the Fourier transform, the so-called short-time Fourier
transform, attempts to circumvent this limitation by applying Fourier
transforms to localized, smoothed time slices of the signal. However, this
leads to limitations in the frequency and time resolution, which the wavelet
transform does not have. Due to its advantages, the wavelet transform is
widely used in the field of image compression (e.g. JPEG 2000). The
superiority of the wavelet transform in this area has accordingly been
demonstrated many times (Santa-Cruz et al., 2002; Woodring et al., 2011; Lewis
& Knowles, 1992). In summary, the wavelet transform allows us to process non-
periodic signals, independent of position and scaling with respect to their
frequency. This allows us to compensate for negative effects that affect
frequency, such as the low-pass effect mentioned above.

To evaluate the effectiveness of our method, we consider challenging super-
resolution problems with complex data sets from physical elasto-plastic
simulations. Super-resolution methods are a common technique to speed up the
necessary calculations and to allow for more control. This can be seen as a
form of post-processing where one simulates only a low-resolution simulation
and uses an up-sampling technique to approximate the behavior of a costly
high-resolution simulation. Neural networks are of special interest here
because of their capability to efficiently approximate the strongly nonlinear
behavior of physical simulations. Applying neural networks to space-time data
sets of physical simulations has seen strongly growing interest in recent
years (Ladicky et al., 2015; Kim et al., 2020), and is particularly
interesting in this context to incorporate additional constraints, e.g., for
temporal coherence (Xie et al., 2018), or for physical plausibility (Tompson
et al., 2017; Kim et al., 2019).

To summarize, the central contributions of our work are: (1) A novel frequency
aware wavelet-based loss formulation, (2) The processing of non-periodic
signals with respect to their frequency, (3) The identification of details
independent of position and scale, (4) A time consistent spatio-temporal up-
sampling of complex physical surfaces. We tested the capabilities of our
approach with three sets of spatio-temporal data from physical simulations.
The experiments show an improvement in synthetic test cases, but the
performance was hardly increased in more complex cases with our setup.

#### Related Work

The idea of exploiting the frequency spectrum of data in Deep Learning is not
new. For example, Mildenhall et al. (2020) and Tancik et al. (2020) use a
Fourier transform as a feature mapping to increase the performance of MLPs.
They show that by providing high-frequency features, the spectral bias can be
reduced, allowing better reconstruction of high-frequency details. The Fourier
transform was also used in the generative domain to reconstruct details from
low-information data (Ichimura, 2018), or to reduce noise in perceptional data
(Yoo et al., 2018; Yadav et al., 2021). Wavelet transformation also found
application in recent publications. For example, wavelets have been used to
augment the provided features to neural networks (Gao & Xiong, 2016; Liu et
al., 2019) as well to variational autoencoders (Ichimura, 2018) and to
generative adversarial networks (Liu et al., 2018; Wang et al., 2020; Gal et
al., 2021). This led to improvements in a variety of generative tasks, such as
denoising (Liu et al., 2020; Zhang et al., 2021), style transfer (Yoo et al.,
2019) or for super-resolution tasks (Zhang et al., 2019). In Gal et al. (2021)
wavelet-based features were provided to the generator and discriminator
architecture, which forces a frequency-aware latent representation. Other
methods attempted to generate wavelet decompositions directly (Huang et al.,
2019; Wang et al., 2020; Zhang et al., 2019).

Most of the generative methods mentioned above focus on perceptual data, such
as photos or videos. However, we focus on structured, noise-free data from
physical simulations, whose properties can be optimally exploited by our
method. Deep learning methods in conjunction with physical models were
employed in variety of contexts, ranging from learning models for physical
intuition (Battaglia et al., 2016; Sanchez-Gonzalez et al., 2018), over
robotic control (Schenck & Fox, 2018; Hu et al., 2019) to engineering
applications (Ling et al., 2016; Morton et al., 2018). In the following, we
focus on fluid-like materials with continuous descriptions, which encompass a
wide range of behavior and pose challenging tasks for learning methods (Mrowca
et al., 2018; Li et al., 2019). For fluid flows in particular, a variety of
learning methods were proposed (Tompson et al., 2017; Prantl et al., 2017; Um
et al., 2018). A common approach to reduce the high computational cost of a
simulation is to employ super-resolution techniques (Dong et al., 2016; Chu &
Thuerey, 2017; Bai et al., 2019). In this context, our work targets the up-
sampling for physics-based animations, similar to the approach proposed by Xie
et al. (2018). However, in contrast to this work, we target phenomena with
clear interfaces, which motivates the frequency-based viewpoint of our work.

For sharp interfaces, Lagrangian models are a very popular discretization of
continuum mechanical systems. E.g., smoothed particle hydrodynamics (SPH)
(Gingold & Monaghan, 1977; Koschier et al., 2019) is a widely-used particle-
based simulation method. While points and particles are likewise frequently
used representations for physical deep learning (Li et al., 2019; Ummenhofer
et al., 2019; Sanchez-Gonzalez et al., 2020), Eulerian, i.e., grid-based
representations offer advantages in terms of efficient and robust kernel
evaluations. A further advantage that results from the observation in wave
space is the ability to process complex time phenomena. The time dimension was
also taken into account in natural imaging works, e.g., by Saito et al. in the
form of a temporal generator (Saito et al., 2017), or via a stochastic
sequence generator (Yu et al., 2017). Other works have included direct $L_{2}$
loss terms as temporal regularizers (Bhattacharjee & Das, 2017; Chen et al.,
2017), which, however, typically strongly restricts the changes over time.
Similar to flow advection, video networks also often use warping information
to align data over time (Liu et al., 2017; de Bezenac et al., 2017). We will
demonstrate that recurrent architectures similar to those used for video
super-resolution (Sajjadi et al., 2018) are likewise very amenable for
physical problems over time.

##  2 Method

The input for our method is a coarsely approximated source simulation, with
the learning objective to infer the surface of a target simulation over space
and time. This target is typically computed via a potentially very costly,
finely resolved simulation run for the same physical setup. When it comes to
the possibilities of simulation representations, there is a great variance. In
our case we have chosen an implicit representation of the data, by a signed-
distance field (SDF) denoted by $g:\mathbb{R}^{3}\rightarrow\mathbb{R}$. An
SDF returns, for a given point, the signed distance to the surface, with
negative being inside the medium. Such a function is realized in practice by a
grid $X\in\mathbb{R}^{M_{x}\times M_{y}\times M_{z}}$, storing the pre-
computed signed distance values, where $M_{*},*\in\\{x,y,z\\}$ specifies the
size of the grid in the respective dimension $x$, $y$ or $z$. We have chosen
this representation because most neural network layers are designed for array-
like representations, and the loss functions on grid-based data are very
efficient to evaluate. Additionally, an implicit representation via a grid can
leverage tools from the field of level-set processing (Adalsteinsson &
Sethian, 1999), and facilitates the frequency viewpoint via a wavelet
transformation. Additional values, like the velocity, are also mapped on a
grid $V\in\mathbb{R}^{M_{x}\times M_{y}\times M_{z}\times 3}$. Our goal is to
let a generative network $\mathcal{G}:\mathbb{R}^{M_{x}\times M_{y}\times
M_{z}\times 4}\rightarrow\mathbb{R}^{N_{x}\times N_{y}\times N_{z}}$ infer a
grid $\tilde{Y}$ which approximates a desired high-resolution simulation
$Y\in\mathbb{R}^{N_{x}\times N_{y}\times N_{z}}$ with $N_{*}=kM_{*},N_{*}\in
N_{\\{x,y,z\\}}$ and up-sample factor $k\in\mathbb{N}$, i.e.
$\mathcal{G}(X)=\tilde{Y}\approx Y$. As our method only requires position and
velocity data from a simulation, it is largely agnostic to the type of solver
or physical model for generating the source and target particle data.

###  2.1 Neural Network Formulation

Figure 1:  The generator takes low-resolution data $X_{t}$ as input, together
with a high-resolution version, from the previous time step $\tilde{Y}_{t-1}$.
The low-resolution data is tri-linearly up-scaled (orange layer), while the
data with the previous state is advected with the velocity of the input (green
layer). The processed inputs are concatenated and processed via four residual
blocks (in blue). A high-resolution frame $\tilde{Y}_{t}$ is generated as
output.

Our method is based on a generative, neural network with a 3D fully-
convolutional ResNet architecture (He et al., 2016) that produces an output
field at a single instance in time. The low-resolution input data is first up-
sampled with a tri-linear up-sampling and then processed with several
convolutional layers, as shown in Figure 1. We use leaky ReLU as activation
function after each layer, except for the last layer, where we use a no
activation. In our case, the input data consists of the implicitly represented
geometry data $X_{t}$, the velocity $V_{t}$ of the simulation as well as the
results of a previous pass $\tilde{Y}_{t}$. The previously generated data is
advected with the low-resolution velocity before further processing. Through
this feedback loop we train our network recurrently by iterating over a
sequence of $T=10$ frames. This yields stability over longer periods of time
and gives better insights about temporal behaviour. Furthermore, the recurrent
training is important to enable persistent behavior over time, such as the
progression of fine surface waves. Unlike the process for generating the input
data, the network training cannot resort to a physical simulation with full
resolution, and hence cannot uniquely determine the evolution of future
states. Therefore, its main learning objective is to capture the dynamics of
the target simulations beyond that basic motion computed with an advection
step. For initialization of the undefined first frame $\tilde{Y}_{-1}$ we use
a tri-linear up-sampled version of the input.

###  2.2 Loss Formulation

Figure 2:  The process of discrete wavelet transform, shown with a chirp
signal. For each level, the given data is processed with a low-pass and a
high-pass filter. The filters are applied along the $x$-axis (to the right),
the $y$-axis (downwards) and the $x$-$y$-axis (diagonally). The result of the
low-pass filter is used to generate the components in the next lower level.
This gives a hierarchical gradation. The results of the high-pass filter are
the wavelet coefficients used for the loss (HF0 - HF2). The chosen
representation clearly shows that lower levels represent low frequencies
(HF2), while higher levels filter high frequencies (HF0).

In the following, we keep the generator network constant to evaluate the
differences for a choice of the loss function to measure the differences
between generated and ground-truth data. The most basic loss function is a
simple mean absolute error (MAE):

$\mathcal{L}_{s}=|Y-\tilde{Y}|.$ (1)

This choice has the central disadvantage that it is ill-suited to measure the
similarity or differences of solutions in a multi-modal settings. Multi-model
here means that multiple solutions for a given input exist. A method that
trains with an MAE loss in this setting will learn the expected value of the
output distribution, i.e, the average of the different solutions. However, the
average is typically not a part of the solution set. Thus, the MAE loss often
does not correspond to the correct distance in solution space, based on
significant factors corresponding to the distribution of the solutions. Our
super-sampling setup is such a problem: Due to the low resolution input, the
high resolution details cannot be determined uniquely, resulting in a variety
of possible solutions when up-sampling. Via physical properties of the
material and its temporal sequence, some solutions can be eliminated, but
nonetheless the space of solutions typically remains infinitely large. If a
MAE loss is used, all such samples from the training data set are simply
averaged to obtain a mean value, so that the result does no longer reflect the
level of detail of the ground-truth data. The MAE loss nevertheless gives a
rough direction, and provides a stable learning target. Hence, we still use it
as a component in the final loss formulation, in combination with a wavelet
loss.

Considering the problem in frequency space, the sought after detail consists
of high-frequency features that cannot be represented by the low-resolution
simulation. The mean-squared distance can be seen in our problem similar to a
low-pass filter. To increase the level of detail again, we have to try to
include higher frequencies in the loss calculation. A simple way would be a
Fourier transformation of the data and a L2 distance in Fourier space. A
Fourier transformation expects periodic data, which can be circumvented with a
short time Fourier transformation. For this one transforms always only a
partial area of the space, smoothed with a window function. Thus one receives
a space-frequency representation, whereby the resolution in the space
dimension depends on the size of the subrange.

A next step is the use of a wavelet transformation. Unlike the Fourier
transform, the wavelet transform also explicitly takes into account the
spatial dimensions. Thus, it gives a larger amount of variance in the spatial
dimensions, which in turn means that spatial details can be isolated and
filtered.

For the wavelet loss we transform data into a space-frequency representation
using a multi-scale discrete wavelet transform
$\mathcal{W}(x)=[w_{i},i\in[0,I)]$:

$\begin{split}HF_{i}&=\langle\psi_{hi},LF_{i-1}\rangle\\\
LF_{i}&=\langle\psi_{li},LF_{i-1}\rangle\\\ LF_{0}&=x,\end{split}$ (2)

where $\psi_{hi}$ and $\psi_{li}$ is a high-pass and low-pass filter of a
wavelet filter bank. In our case, we use Daubechies-2 wavelets (Daubechies,
1992). The first pass starts processing the input data $x$. The high pass
filter generates the wavelet component of the respective level, the low pass
filter on the other side generates a scaled down version of the data, which is
then used in the next pass. The process is repeated to generate the wavelet
components in a hierarchical fashion. Each level represents a different
frequency, resulting in a list of differently scaled wavelet components.
Figure 2 shows an example of the discrete wavelet transformation of a 2D cirp
signal. A signal whose frequency increases progressively over time. The number
of levels $I$ is given by $\lfloor log_{2}d\rfloor$, where $d$ is the smallest
input dimension. In addition we scale the values using the logarithm:

$w_{i}=log_{2}(|HF_{i}|+\epsilon).$ (3)

This way we compensate for the fact that high-frequency components have less
energy, making low-frequency components very dominant. The added $\epsilon$ is
used to avoid very high values and to keep the scaling stable. For multi-
dimensional grids, the wavelet transform is applied multiple times, with each
pass applying the filter to the respective dimension. We denote the wavelet
transform by $\mathcal{W}_{d}$ where $d$ denote the dimensions from which the
wavelet transform is formed. Using the wavelet transform, the L1 distance is
calculated from the scaled and transformed data and used as the loss:

$L_{ws}=\sum\limits_{t=0}^{T}|\mathcal{W}_{xyz}(Y_{t})-\mathcal{W}_{xyz}(\tilde{Y}_{t})|$
(4)

While we have primarily focused on spatial content so far, i.e., the surface
of the material, the temporal behavior likewise plays a crucial role, and
poses similar difficulties in our multi-modal setting. On the one hand, the
generation of details can quickly lead to temporally incoherent results, which
is characterized by unappealing flickering. On the other hand, our network
also should be able to match and recreate spatial solutions over time that
reflect the physical behavior. Similar to the spatial loss, we can also resort
to a wavelet transformation. One possible solution would be, to increase the
dimensionality of the wavelet transform by one, giving a space-time-frequency
representation. However, we chose to consider space and time separately. The
separation of spatial and temporal frequencies gives us more flexibility by
adjusting the weighting. Using the wavelet transform, we generate a time-
frequency representation, similar to the spatial loss:

$L_{wt}=|\mathcal{W}_{t}(Y)-\mathcal{W}_{t}(\tilde{Y})|.$ (5)

The distances in wavelet space are included in the loss formulation of the
generator which gives the final loss function:

$\mathcal{L}_{G}=\mathcal{L}_{s}+\alpha\mathcal{L}_{ws}+\beta\mathcal{L}_{wt},$
(6)

where $\alpha$ and $\beta$ indicate the weighting of the individual loss
terms.

To indicate the focus on surface structures, we refer to the final version of
our generative network as surfNet. For a more detailed description of the
training and the network architecture we refer to the appendix A.1 and A.2.

\begin{overpic}[width=212.47617pt,trim=0.0pt 20.075pt 0.0pt
50.18748pt]{data/spatial_0029.png}
\put(-8.0,15.0){\rotatebox{90.0}{\color[rgb]{0,0,0}{amplitude}}}
\put(30.0,-10.0){\color[rgb]{0,0,0}{frequency [Hz]}}\end{overpic}

\begin{overpic}[width=212.47617pt,trim=0.0pt 20.075pt 0.0pt
50.18748pt]{data/temporal_0029.png}
\put(-8.0,15.0){\rotatebox{90.0}{\color[rgb]{0,0,0}{amplitude}}}
\put(30.0,-10.0){\color[rgb]{0,0,0}{frequency [Hz]}}\end{overpic}

(a) Frequency histograms (spatial and temp.).

(b) Synthetic data sample.

Figure 3:  The averaged frequency spectrum of our data set (a) for the
surface, spatial (left), and for the temporal behavior (right). The spectrum
of the low frequency input data is shown in blue, while in orange the spectrum
of the ground-truth data is shown. In the spatial spectrum, it’s clearly
visible that the spectrum of the input data only covers a quarter of the
frequencies, while a large peak of the ground-truth is visible for higher
frequencies. This highlights what the generator needs to reproduce. A similar
shape can be seen in the time spectrum as the time response is strongly
coupled to the surface frequency. Here, the spectrum of the input data (blue)
continues as the time discretization is the same for both. Figure (b) shows an
example of a surface from our synthetic data set: the input wave in blue, and
the corresponding ground-truth target in purple.

##  3 Results

For evaluation we consider three different data sets: a synthetic 2D case, a
2D simulation based on a mass-spring system and a 3D particle-based
simulation. The synthetic data is fully controllable such that the frequency
spectrum of the surface can be evaluated reliably. The basis is a wavy
surface, based on a sine wave with varying frequency (Figure 3). This forms a
wide range of analysis to isolate problems in the generative process and
illustrate the aspects of the proposed method. We then evaluate the
established methodology with data generated from physical simulations for a
more complex scenario. We use elastic bodies based on a mass-spring system in
2D and we also consider a simulation with plastic material with data generated
from a highly viscous SPH simulation (Weiler et al., 2018). For more
information on the three data sets, refer to the appendix A.3, A.4 and A.5.

\begin{overpic}[width=397.48499pt,trim=50.18748pt 50.18748pt 50.18748pt
50.18748pt,clip]{figures/freq/baseline_spatial.pdf}
\put(-8.0,20.0){\rotatebox{90.0}{\color[rgb]{0,0,0}{log.
amplitude}}}\end{overpic}\begin{overpic}[width=397.48499pt,trim=50.18748pt
50.18748pt 50.18748pt 50.18748pt,clip]{figures/freq/baseline_temporal.pdf}
\put(-8.0,20.0){\rotatebox{90.0}{\color[rgb]{0,0,0}{log. amplitude}}}
\put(30.0,-8.0){\color[rgb]{0,0,0}{frequency [Hz]}}\end{overpic}

(a) MAE.

\begin{overpic}[width=397.48499pt,trim=50.18748pt 50.18748pt 50.18748pt
50.18748pt,clip]{figures/freq/rfft_temporal.pdf}
\put(30.0,-8.0){\color[rgb]{0,0,0}{frequency [Hz]}}\end{overpic}

(b) RFFT loss.

\begin{overpic}[width=397.48499pt,trim=50.18748pt 50.18748pt 50.18748pt
50.18748pt,clip]{figures/freq/gan_temporal.pdf}
\put(30.0,-8.0){\color[rgb]{0,0,0}{frequency [Hz]}}\end{overpic}

(c) GAN.

\begin{overpic}[width=397.48499pt,trim=50.18748pt 50.18748pt 50.18748pt
50.18748pt,clip]{figures/freq/wvlt_temporal.pdf}
\put(30.0,-8.0){\color[rgb]{0,0,0}{frequency [Hz]}}\end{overpic}

(d) surfNet.

Figure 4:  Frequency spectrum comparison of different versions of our network for the synthetic 2D data set. The top row always shows the spatial frequency, while the bottom row shows the temporal frequency. A logarithmic scaling is used for the amplitude. The spectrum of the input is shown in blue, the ground-truth in pink, and the predictions in green. Version (a) was trained with a simple MAE loss, while (b) was extended with a RFFT loss; (c) was trained with an adversarial loss and (d) was trained with our wavelet loss.  |  MAE | RFFT | GAN |  surfNet
---|---|---|---|---
MAE | 0.039 | 0.074 | 0.052 | 0.047
Spatial freq. MAE | 0.990 | 1.110 | 0.819 | 0.884
Temp. freq. MAE | 0.881 | 0.997 | 0.978 | 0.738
Table 1:  Mean errors for 2D variants. The spatial and temporal frequency MAE
values represent the averaged difference of the generated and the target
residuals in frequency space.

\begin{overpic}[width=433.62pt,trim=642.4pt 200.74998pt 642.4pt
401.49998pt,clip]{data/2d/wlt/110.jpg}
\put(-50.0,50.0){\color[rgb]{0,0,0}{$a)$}}\end{overpic}

\begin{overpic}[width=433.62pt,trim=1264.72499pt 200.74998pt 20.075pt
401.49998pt,clip]{data/2d/mse/110.jpg}
\put(-50.0,50.0){\color[rgb]{0,0,0}{$b)$}}\end{overpic}

\begin{overpic}[width=433.62pt,trim=1264.72499pt 200.74998pt 20.075pt
401.49998pt,clip]{data/2d/wlt/110.jpg}
\put(-50.0,50.0){\color[rgb]{0,0,0}{$c)$}}\end{overpic}

(a)

(b) Spatial frequency.

\begin{overpic}[width=433.62pt,trim=642.4pt 200.74998pt 642.4pt
401.49998pt,clip]{data/2d/wlt/110.jpg}
\put(-50.0,50.0){\color[rgb]{0,0,0}{$a)$}}\end{overpic}

\begin{overpic}[width=433.62pt,trim=1264.72499pt 200.74998pt 20.075pt
401.49998pt,clip]{data/2d/mse/110.jpg}
\put(-50.0,50.0){\color[rgb]{0,0,0}{$b)$}}\end{overpic}

\begin{overpic}[width=433.62pt,trim=1264.72499pt 200.74998pt 20.075pt
401.49998pt,clip]{data/2d/wlt/110.jpg}
\put(-50.0,50.0){\color[rgb]{0,0,0}{$c)$}}\end{overpic}

(c) Data.

(d) Temporal frequency.

Figure 5:  A comparison of wavelet-based features between ground-truth data
(a) and the generated data with MAE loss (b) and our surfNet (c). For this,
the spatial spectrum of a selected frame was generated, together with the
temporal spectrum for the corresponding sequence. The respective columns
implicitly indicate the frequency (from left, higher to right, lower
frequency). This clearly shows how the high-frequency features are revealed by
the wavelet transform and partially reconstructed by our method.

\begin{overpic}[width=433.62pt,trim=50.18748pt 281.04999pt 50.18748pt
451.68748pt,clip]{data/2d/wlt/64.jpg}
\put(12.0,1.5){\color[rgb]{1,1,1}{{Source}}}
\put(43.0,1.5){\color[rgb]{1,1,1}{{Ground-truth}}}
\put(80.0,1.5){\color[rgb]{1,1,1}{{surfNet}}}\end{overpic}

\begin{overpic}[width=433.62pt,trim=1264.72499pt 281.04999pt 50.18748pt
451.68748pt,clip]{data/2d/mse/64.jpg}
\put(45.0,5.0){\color[rgb]{1,1,1}{{MAE}}}\end{overpic}

\begin{overpic}[width=433.62pt,trim=1264.72499pt 281.04999pt 50.18748pt
451.68748pt,clip]{data/2d/freq/64.jpg}
\put(40.0,5.0){\color[rgb]{1,1,1}{{RFFT}}}\end{overpic}

\begin{overpic}[width=433.62pt,trim=1264.72499pt 281.04999pt 50.18748pt
451.68748pt,clip]{data/2d/gan/64.jpg}
\put(42.0,5.0){\color[rgb]{1,1,1}{{GAN}}}\end{overpic}

Figure 6:  Comparison of a frame from a test run with the synthetic 2D data
set. The input is shown in blue, ground-truth in red, and predictions in
green.

###  3.1 Surface Waves

We first consider the controlled 2D data. Training is performed with the full
data set of 10000 samples with 10 frames, and we evaluate the resulting models
on 10 simulations that are not part of the training set with 40 frames each.
The data consists of one-dimensional sine waves, from which a 2D surface was
generated. This intuitive setup allows us to easily evaluate the spectrum of
the surface.

In the following evaluations we compare our approach with a Fourier and a GAN-
based loss formulation. The Fourier loss works similarly to our wavelet Loss,
except that we use a real Fourier transform (RFFT) instead of the wavelet
transform. The RFFT version serves mainly as a sanity check, since the setup
is based on a combination of periodic sine waves and thus falls exactly within
the scope of the Fourier transform. This makes the task much easier for the
RFFT. With more complex problems, the Fourier transformation comes to its
limits, since the signals no longer occur periodically. Thus, our goal is to
achieve comparable results with the more flexible wavelet loss. The comparison
with the GAN version (Xie et al., 2018), however, serves to show that our
method is similarly powerful while retaining the advantage that the loss
function is interpretable. As part of our evaluation, we analyse the frequency
of the surface. In Figure 4, we compare the spectrum of the data generated by
the networks with those from the ground truth. Figure 4(a) shows the spectrum
of generated data of a variant of our network trained with a MAE loss. It
hardly shows any deflection in the high frequency range. If we now use our
wavelet loss (Figure 4(d)), we can clearly see that higher frequencies are
more present in the generated data, making the spectrum similar to that of the
ground-truth. In Figure 4(b) and Figure 4(c) you can see the spectrum obtained
with the RFFT setup and the GAN setup.

In Table 1 we also compare the mean absolute error values of the data
quantitatively in Euclidean space and in frequency space. The mean errors
likewise illustrate that the MAE version is not able to reconstruct the high
frequencies of the target function.

Another illustration of the superiority of our method over the baseline is
shown in Figure 5. It shows the wavelet transformation from a chosen sample.
From top to bottom, it compares ground-truth, a network with MAE loss, and the
surfNet result. While the MAE-based variant is not able to reconstruct fine
details, the surfNet can recover these, as highlighted by the wavelet
representation. The high-frequency blocks clearly illustrate the differences
into terms of reconstruction accuracy.

As a qualitative evaluation, Figure 6 shows visual examples of the 2D test
data set. While a version based on MAE generates very smooth images in the
random data set, the surfNet versions are able to reconstruct the jagged
edges. In direct comparison with the ground-truth data there are differences,
but this is due to the randomness of the data. Therefore, the exact solution
cannot be reconstructed, but the surfNet is still able to generate a very
plausible solution.

###  3.2 Mass-Spring System

| MAE | surfNet
---|---|---
MAE | 0.017 | 0.020
Spatial freq. MAE | 0.813 | 0.809
Temp. freq. MAE | 0.523 | 0.466
Table 2:  Mean error measurements for 2D test scenarios.

For a more complex test, we considered the behavior of different elastic 2D
shapes simulated on the basis of mass-spring systems. The low-resolution
version, which serves as input, provides the rough motion with a static shape.
The network has the task to reconstruct the oscillations caused by the elastic
material. As you can see in Figure 7, the generator is able to reconstruct the
details to some extent for both variants. The results with the wavelet loss
(Figure 7) have, as expected, more high-frequency features, while they are
more subtle with MAE. Nevertheless, the generated details of our method do not
always seem to match the physical behavior of the data. Apart from that, the
differences to the ground truth are still very large. It seems like the
network in this case is overfitting to certain high-frequency details rather
than generalizing to the underlying physical behavior.

(a) Frame 45.

(a) Frame 66.

(a) Frame 123 surfNet.

(b) Frame 123 MAE.

Figure 7:  Comparison of 3D frames with a surfNet model (a) and with a MAE-
based generator (b). The input is shown in blue, ground-truth in red, and
predictions in green.

(a) Frame 98.

(a) Frame 128.

(a) Frame 164 surfNet.

(b) Frame 164 MAE.

Figure 8:  Comparison of 3D frames with a surfNet model (a) and with a MAE-
based generator (b). The input is shown in blue, ground-truth in red, and
predictions in green.

###  3.3 SPH Simulation

| MAE | surfNet
---|---|---
MAE | 0.087 | 0.107
Spatial freq. MAE | 0.962 | 0.695
Temp. freq. MAE | 0.464 | 0.250
Table 3:  Mean error measurements for 3D test scenarios.

A further test case based on physical simulation we use data based on a high
viscosity 3D SPH simulation. We compare our method with a MAE based approach,
whose is not able to reconstruct the outlines and details of the ground-truth
data and is temporally incoherent. Our method achieves slightly better
results, as can also be seen in Table 3. Both, errors in the spatial and
temporal frequencies are smaller for our method.

Figure 8(a) shows that our prediction is able to reconstruct and preserve most
of the details, even if they are not present in the input. We have
deliberately chosen frames after a long run-time (98-164 time-steps) to show
that details can persist and that they are the result of complex, physically-
based behavior over multiple frames. Again, we compared our method with a
simpler MAE-based generator. With MAE, the results tend to have unsightly
artifacts that lead to flickering over time, as can be seen in Figure 8(b).
The wavelet variant (Figure 8(c)) delivers much more coherent results, but
similar to the 2D case, the amount of additional detail was not much higher.

##  4 Conclusion

Our method sets the foundation for a generic, frequency-based loss function
that can find application in many spatial and temporal geometric learning
problems. This is interesting especially for physical problems due to the
structural nature of the data. The effectiveness of our method has been well
demonstrated by evaluations on a synthetic data set. The evaluation of more
complex problems was less successful, but we are convinced that due to the
generic and simple nature of our method it can still be useful for more
complex problems if applied correctly. Our method provides a first step
towards evaluation and synthesis of physical space-time processes, and could
be employed for other phenomena such as turbulence Ling et al. (2016) or
weather Zaytar & El Amrani (2016). Furthermore, it will be interesting to
employ it in conjunction with other frequency-based representations Sitzmann
et al. (2020).

## References

  * Adalsteinsson & Sethian (1999) David Adalsteinsson and James A Sethian.  The fast construction of extension velocities in level set methods.  _Journal of Computational Physics_ , 148(1):2–22, 1999.
  * Bai et al. (2019) Kai Bai, Wei Li, Mathieu Desbrun, and Xiaopei Liu.  Dynamic upsampling of smoke through dictionary-based learning.  _arXiv:1910.09166_ , 2019.
  * Battaglia et al. (2016) Peter W. Battaglia, Razvan Pascanu, Matthew Lai, Danilo Rezende, and Koray Kavukcuoglu.  Interaction networks for learning about objects, relations and physics.  In _Advances in Neural Information Processing Systems_ , 2016.
  * Bender (2017) Jan Bender.  Splishsplash library, 2017.
  * Bhattacharjee & Das (2017) Prateep Bhattacharjee and Sukhendu Das.  Temporal coherency based criteria for predicting video frames using deep multi-stage generative adversarial networks.  In _Advances in Neural Information Processing Systems_ , pp. 4271–4280, 2017.
  * Bruna et al. (2016) Joan Bruna, Pablo Sprechmann, and Yann LeCun.  Super-resolution with deep convolutional sufficient statistics.  In Yoshua Bengio and Yann LeCun (eds.), _4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings_ , 2016.  URL http://arxiv.org/abs/1511.05666.
  * Chen et al. (2017) Dongdong Chen, Jing Liao, Lu Yuan, Nenghai Yu, and Gang Hua.  Coherent online video style transfer.  In _The IEEE International Conference on Computer Vision (ICCV)_ , Oct 2017.
  * Chu & Thuerey (2017) Mengyu Chu and Nils Thuerey.  Data-driven synthesis of smoke flows with CNN-based feature descriptors.  _ACM Trans. Graph._ , 36(4)(69), 2017.
  * Daubechies (1992) Ingrid Daubechies.  _Ten lectures on wavelets_.  SIAM, 1992.
  * de Bezenac et al. (2017) Emmanuel de Bezenac, Arthur Pajot, and Patrick Gallinari.  Deep learning for physical processes: Incorporating prior scientific knowledge.  _arXiv preprint arXiv:1711.07970_ , 2017.
  * Dong et al. (2016) Chao Dong, Chen Change Loy, Kaiming He, and Xiaoou Tang.  Image super-resolution using deep convolutional networks.  _IEEE transactions on pattern analysis and machine intelligence_ , 38(2):295–307, 2016.
  * Dosovitskiy & Brox (2016) Alexey Dosovitskiy and Thomas Brox.  Generating images with perceptual similarity metrics based on deep networks.  _CoRR_ , abs/1602.02644, 2016.  URL http://arxiv.org/abs/1602.02644.
  * Gal et al. (2021) Rinon Gal, Dana Cohen, Amit Bermano, and Daniel Cohen-Or.  SWAGAN: A style-based wavelet-driven generative model.  _CoRR_ , abs/2102.06108, 2021.  URL https://arxiv.org/abs/2102.06108.
  * Gao & Xiong (2016) Xing Gao and Hongkai Xiong.  A hybrid wavelet convolution network with sparse-coding for image super-resolution.  In _2016 IEEE International Conference on Image Processing (ICIP)_ , pp. 1439–1443, 2016.  doi: 10.1109/ICIP.2016.7532596.
  * Gingold & Monaghan (1977) Robert A Gingold and Joseph J Monaghan.  Smoothed particle hydrodynamics: Theory and application to non-spherical stars.  _Monthly Notices of the Royal Astronomical Society_ , 181(3), 1977.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proc. of IEEE Comp. Vision and Pattern Rec._ , pp. 770–778, 2016.
  * Hu et al. (2019) Yuanming Hu, Luke Anderson, Tzu-Mao Li, Qi Sun, Nathan Carr, Jonathan Ragan-Kelley, and Frédo Durand.  Difftaichi: Differentiable programming for physical simulation.  _arXiv preprint arXiv:1910.00935_ , 2019.
  * Huang et al. (2019) Huaibo Huang, Ran He, Zhenan Sun, and Tieniu Tan.  Wavelet domain generative adversarial network for multi-scale face hallucination.  _International Journal of Computer Vision_ , 127:763–784, 2019.
  * Ichimura (2018) Naoyuki Ichimura.  Spatial frequency loss for learning convolutional autoencoders.  _CoRR_ , abs/1806.02336, 2018.  URL http://arxiv.org/abs/1806.02336.
  * Johnson et al. (2016) Justin Johnson, Alexandre Alahi, and Li Fei-Fei.  Perceptual losses for real-time style transfer and super-resolution.  In _European Conference on Computer Vision_ , pp. 694–711. Springer, 2016.
  * Kim et al. (2019) Byungsoo Kim, Vinicius C. Azevedo, Nils Thuerey, Theodore Kim, Markus Gross, and Barbara Solenthaler.  Deep fluids: A generative network for parameterized fluid simulations.  _Computer Graphics Forum_ , 2019.  ISSN 1467-8659.  doi: 10.1111/cgf.13619.
  * Kim et al. (2020) Byungsoo Kim, Vinicius C Azevedo, Markus Gross, and Barbara Solenthaler.  Lagrangian neural style transfer for fluids.  _ACM Trans. Graph._ , 39(4), 2020.
  * Koschier et al. (2019) Dan Koschier, Jan Bender, Barbara Solenthaler, and Matthias Teschner.  Smoothed particle hydrodynamics techniques for the physics based simulation of fluids and solids.  In _EUROGRAPHICS 2019 Tutorials_. Eurographics Association, 2019\.
  * Ladicky et al. (2015) Lubor Ladicky, SoHyeon Jeong, Barbara Solenthaler, Marc Pollefeys, and Markus Gross.  Data-driven fluid simulations using regression forests.  _ACM Trans. Graph._ , 34(6):199, 2015.
  * Lewis & Knowles (1992) A.S. Lewis and G. Knowles.  Image compression using the 2-d wavelet transform.  _IEEE Transactions on Image Processing_ , 1(2):244–250, 1992.  doi: 10.1109/83.136601.
  * Li et al. (2019) Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua B. Tenenbaum, and Antonio Torralba.  Learning particle dynamics for manipulating rigid bodies, deformable objects, and fluids.  In _International Conference on Learning Representations_ , 2019.
  * Ling et al. (2016) Julia Ling, Andrew Kurzawski, and Jeremy Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  _Journal of Fluid Mechanics_ , 807, 2016.
  * Liu et al. (2017) Ding Liu, Zhaowen Wang, Yuchen Fan, Xianming Liu, Zhangyang Wang, Shiyu Chang, and Thomas Huang.  Robust video super-resolution with learned temporal dynamics.  In _The IEEE International Conference on Computer Vision (ICCV)_ , Oct 2017.
  * Liu et al. (2020) Lin Liu, Jianzhuang Liu, Shanxin Yuan, Gregory G. Slabaugh, Ales Leonardis, Wengang Zhou, and Qi Tian.  Wavelet-based dual-branch network for image demoireing.  _CoRR_ , abs/2007.07173, 2020.  URL https://arxiv.org/abs/2007.07173.
  * Liu et al. (2019) Pengju Liu, Hongzhi Zhang, Wei Lian, and Wangmeng Zuo.  Multi-level wavelet convolutional neural networks.  _IEEE Access_ , 7:74973–74985, 2019.  doi: 10.1109/ACCESS.2019.2921451.
  * Liu et al. (2018) Yunfan Liu, Qi Li, and Zhenan Sun.  Attribute enhanced face aging with wavelet-based generative adversarial networks.  _CoRR_ , abs/1809.06647, 2018.  URL http://arxiv.org/abs/1809.06647.
  * Mildenhall et al. (2020) Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik, Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  _CoRR_ , abs/2003.08934, 2020.  URL https://arxiv.org/abs/2003.08934.
  * Morton et al. (2018) Jeremy Morton, Antony Jameson, Mykel J Kochenderfer, and Freddie Witherden.  Deep dynamical modeling and control of unsteady fluid flows.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Mrowca et al. (2018) Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haber, Li F Fei-Fei, Josh Tenenbaum, and Daniel L Yamins.  Flexible neural representation for physics prediction.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Prantl et al. (2017) Lukas Prantl, Boris Bonev, and Nils Thuerey.  Pre-computed liquid spaces with generative neural networks and optical flow.  _arXiv:1704.07854_ , 2017.
  * Saito et al. (2017) Masaki Saito, Eiichi Matsumoto, and Shunta Saito.  Temporal generative adversarial nets with singular value clipping.  In _IEEE International Conference on Computer Vision (ICCV)_ , pp. 2830–2839, 2017.
  * Sajjadi et al. (2018) Mehdi SM Sajjadi, Raviteja Vemulapalli, and Matthew Brown.  Frame-recurrent video super-resolution.  In _IEEE Conference on Computer Vision and Pattern Recognition_ , 2018\.
  * Sanchez-Gonzalez et al. (2018) Alvaro Sanchez-Gonzalez, Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin A. Riedmiller, Raia Hadsell, and Peter W. Battaglia.  Graph networks as learnable physics engines for inference and control.  In _International Conference on Machine Learning_ , 2018.
  * Sanchez-Gonzalez et al. (2020) Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter W Battaglia.  Learning to simulate complex physics with graph networks.  In _International Conference on Machine Learning_ , 2020.
  * Santa-Cruz et al. (2002) Diego Santa-Cruz, Raphaël Grosbois, and Touradj Ebrahimi.  Jpeg 2000 performance evaluation and assessment.  _Signal Processing: Image Communication_ , 17(1):113–130, 2002.  ISSN 0923-5965.  doi: https://doi.org/10.1016/S0923-5965(01)00025-X.  URL https://www.sciencedirect.com/science/article/pii/S092359650100025X.  JPEG 2000.
  * Schenck & Fox (2018) Connor Schenck and Dieter Fox.  SPNets: Differentiable fluid dynamics for deep neural networks.  In _Conference on Robot Learning_ , 2018.
  * Sitzmann et al. (2020) Vincent Sitzmann, Julien NP Martel, Alexander W Bergman, David B Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  _arXiv preprint arXiv:2006.09661_ , 2020.
  * Tancik et al. (2020) Matthew Tancik, Pratul P. Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T. Barron, and Ren Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  _CoRR_ , abs/2006.10739, 2020.  URL https://arxiv.org/abs/2006.10739.
  * Tompson et al. (2017) Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin.  Accelerating eulerian fluid simulation with convolutional networks.  _ICML_ , 2017.
  * Um et al. (2018) Kiwon Um, Xiangyu Hu, and Nils Thuerey.  Liquid splash modeling with neural networks.  In _Computer Graphics Forum_ , volume 37(8), pp. 171–182. Wiley Online Library, 2018.
  * Ummenhofer et al. (2019) Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun.  Lagrangian fluid simulation with continuous convolutions.  In _International Conference on Learning Representations_ , 2019.
  * Wang et al. (2020) Jianyi Wang, Xin Deng, Mai Xu, Congyong Chen, and Yuhang Song.  Multi-level wavelet-based generative adversarial network for perceptual quality enhancement of compressed video.  In Andrea Vedaldi, Horst Bischof, Thomas Brox, and Jan-Michael Frahm (eds.), _Computer Vision – ECCV 2020_ , pp. 405–421, Cham, 2020. Springer International Publishing.  ISBN 978-3-030-58568-6.
  * Weiler et al. (2018) Marcel Weiler, Dan Koschier, Magnus Brand, and Jan Bender.  A physically consistent implicit viscosity solver for sph fluids.  _Computer Graphics Forum_ , 37(2), 2018.
  * Woodring et al. (2011) Jonathan Woodring, Susan Mniszewski, Christopher Brislawn, David DeMarle, and James Ahrens.  Revisiting wavelet compression for large-scale climate data using jpeg 2000 and ensuring data precision.  In _2011 IEEE Symposium on Large Data Analysis and Visualization_ , pp. 31–38, 2011.  doi: 10.1109/LDAV.2011.6092314.
  * Xie et al. (2018) You Xie, Erik Franz, Mengyu Chu, and Nils Thuerey.  tempogan: A temporally coherent, volumetric gan for super-resolution fluid flow.  _ACM Trans. Graph._ , 37(4):95:1–95:15, 2018\.  doi: 10.1145/3197517.3201304.
  * Yadav et al. (2021) Ojasvi Yadav, Koustav Ghosal, Sebastian Lutz, and Aljosa Smolic.  Frequency-domain loss function for deep exposure correction of dark images.  _Signal, Image and Video Processing_ , May 2021.  ISSN 1863-1711.  doi: 10.1007/s11760-021-01915-4.  URL http://dx.doi.org/10.1007/s11760-021-01915-4.
  * Yoo et al. (2019) Jaejun Yoo, Youngjung Uh, Sanghyuk Chun, Byeongkyu Kang, and Jung-Woo Ha.  Photorealistic style transfer via wavelet transforms.  _CoRR_ , abs/1903.09760, 2019.  URL http://arxiv.org/abs/1903.09760.
  * Yoo et al. (2018) Jaeyoung Yoo, Sang-Ho Lee, and Nojun Kwak.  Image restoration by estimating frequency distribution of local patches.  _CoRR_ , abs/1805.09097, 2018.  URL http://arxiv.org/abs/1805.09097.
  * Yu et al. (2017) Lantao Yu, Weinan Zhang, Jun Wang, and Yong Yu.  Seqgan: Sequence generative adversarial nets with policy gradient.  In _AAAI_ , pp. 2852–2858, 2017.
  * Zaytar & El Amrani (2016) Mohamed Akram Zaytar and Chaker El Amrani.  Sequence to sequence weather forecasting with long short-term memory recurrent neural networks.  _International Journal of Computer Applications_ , 143(11):7–11, 2016.
  * Zhang et al. (2019) Qi Zhang, Huafeng Wang, Tao Du, Sichen Yang, Yuehai Wang, Zhiqiang Xing, Wenle Bai, and Yang Yi.  Super-resolution reconstruction algorithms based on fusion of deep learning mechanism and wavelet.  In _Proceedings of the 2nd International Conference on Artificial Intelligence and Pattern Recognition_ , AIPR ’19, pp. 102–107, New York, NY, USA, 2019. Association for Computing Machinery.  ISBN 9781450372299.  doi: 10.1145/3357254.3358600.  URL https://doi.org/10.1145/3357254.3358600.
  * Zhang et al. (2021) Ran Zhang, Junchi Bin, Zheng Liu, and Erik Blasch.  Chapter 13 - wggan: A wavelet-guided generative adversarial network for thermal image translation.  In Arun Solanki, Anand Nayyar, and Mohd Naved (eds.), _Generative Adversarial Networks for Image-to-Image Translation_ , pp. 313–327. Academic Press, 2021.  ISBN 978-0-12-823519-5.  doi: https://doi.org/10.1016/B978-0-12-823519-5.00015-4.  URL https://www.sciencedirect.com/science/article/pii/B9780128235195000154.

##  Appendix A Appendix

###  A.1 Implementation Details

The main part of the generator consists of four blocks containing two
convolution layers and a residual connection. The feature count of the
convolution layers per block is as follows: $[8,16],[32,32],[16,8],[1,1]$. The
kernel size is $5\times 5\times 5$ for all layers. It is important to mention
that for the first convolution of our network we do not use zero padding as
usual but mirror padding. This is because we work with tiles from the training
data and not with the complete frame and zero padding falsifies the values at
the transitions between the tiles.

###  A.2 Training Details

For the training we have implemented our network with the Tensorflow
Framework. We use an Adam optimizer with a learning rate from 0.0005 to
0.00001 with piece-wise constant decay and a batch size of 16 for the 2D tests
and 4 for the 3D tests for 50k iterations. All other weights are initialized
with the respective standard initializers of Tensorflow version 2.1. The
weighting factors $\alpha$ and $\beta$ of Equation 6 are set to $100.0$ and
$10.0$ correspondingly.

Figure 9:  In blue a wave with a randomly varying low frequency $f_{l}$, which
can still be represented by the low-resolution sampling and serves as input
data set for our synthetic data set. In red, on the other hand, we see a wave
with a high frequency $f_{h}$ that can only be correctly represented with much
higher resolution. Combined with the low-frequency version, we get the ground-
truth data (violet) for the synthetic data set.

###  A.3 Synthetic Data Set

We use a synthetic data set to test and evaluate different aspects of our
approach. The data set is designed to have a simple, clearly defined behavior,
and such that the frequency spectrum of the surface can be evaluated reliably.
Therefore, we use a horizontal wavy surface based on a 1D sine function
$s_{t}(x)$ with randomly varying frequency $f_{l}\in(0,\frac{M}{2})$:

$s_{0}(x)=sin(\frac{\pi f_{l}(x)x}{M}),x\in[0,M),$ (7)

where $M$ is the resolution of the low-resolution data.

To make a time sequence out of this, we use a simple wave equation:

$\frac{\delta^{2}s}{\delta t^{2}}=\frac{\delta^{2}s}{\delta x^{2}},$ (8)

Discretizing this we can calculate the vertical velocity of our surface as
follows:

$v_{t}(x)=v_{t-1}(x)+\frac{2*s_{t}(x)-s_{t}(x-\Delta x)+s_{t}(x+\Delta
x)}{\Delta x},$ (9)

where $v_{0}(x)=0$. Given the velocity we can now calculate the next frame as
follows:

$s_{t}(x)=s_{t-\Delta t}(x)+\Delta tv_{t}(x),$ (10)

For the high resolution data set, i.e, the targets to be learned, we use the
same low resolution wave as base, and modulate it with a high frequency
component (Figure 9):

$t_{0}(x)=sin(\frac{\pi f_{l}(x)x}{kM})+sin(\frac{\pi
f_{h}(x)x}{kM}),x\in[0,kM),$ (11)

where $k$ is the chosen up-sampling factor and the frequency $f_{h}(x)$ is
chosen so that it cannot be represented by the low-resolution version.
According to the Nyquist-Shannon sampling theorem, the frequency should be
higher than $\frac{M}{2}$ and below $\frac{kM}{2}$ (Figure 3(a)). The
generation of a sequence is done in the same way like for the low-resolution
data.

Based on this setup, we generate two different data sets of high resolution
data: one where we modulate with a fixed high frequency $f_{h}(x)=const$,
whereas in the second version we vary the this high frequency component
(Figure 3(b)). Thus the first version represents a deterministic up-sampling,
which a generator should be able to reconstruct perfectly, whereas the second
version is ill-posed, i.e. several solutions are possible, and hence poses a
much more difficult learning target. Both data sets consist of 10000 sequences
with 10 frames each at the end. While the deterministic version only serves as
a sanity check, the second, randomized version shows how well the method can
approximate the ground-truth distribution for ill-posed tasks like the actual
super-resolution problem for physical simulation data.

###  A.4 Mass-Spring System Data

As a first, complex, physics-based example, we use data from an 2D elasto-
plastic simulation with a mass-spring system. For this, we generate mass
particles in star-shaped shapes with different scale and rotation and connect
them with springs. We then simulate how the generated objects fall to the
ground. We map the particle-based simulations onto a grid to make them
compatible with the network. We use a high-resolution grid for our ground-
truth data and a low-resolution grid for our input data. The network thus has
to learn to reconstruct from the low-resolution data the high-frequency
details generated at the surface of the ground-truth simulation. For the data,
900 particles are simulated over 400 time steps in 1000 different setups. The
generated grids have a size of $64\times 64$ for the high resolution and
$32\times 32$ for the low resolution.

###  A.5 Simulation Data

For the generation of simulation data we use an SPH solver of the SPlisHSPlasH
framework (Bender, 2017). There are different materials to choose from. With
materials that exhibit high-frequency physical behavior, chaotic behavior
occurs in some cases, such as splashes in water, which are typically very
difficult to reconstruct. With more viscous materials, such as gel, details
are mainly distinguished by folds and fine waves on the surface. The chaotic
behaviour is very difficult to reconstruct and it is sometimes very difficult
to understand how correct the behaviour is. For these reasons we focus more on
materials like gel. With gel, fine wrinkles can form on the surface which
allows a good evaluation of the method. Another special feature is that
details are persistent over time. Methods such as Xie et al. (2018) cannot
represent such details because they do not provide feedback in the network.
This means there is no memory. Therefore, we use an plastic material with high
viscosity, simulated with an advanced viscosity solver (Weiler et al., 2018).

For the training we generate 60 different scenes with 300 frames per
simulation. The time step corresponds to 100 frames per second. The scenes
consist of randomly generated shapes that fall into a pool from different
heights at random times. This creates interesting waves and folds on the
surface. The ground-truth resolution is $160^{3}$. For the generation of the
low-resolution data the ground-truth data is scaled down by the desired up-
sampling factor $k$ and then smoothed with a Gaussian blur. This results in
synchronous data pairs that can be used in a supervised setup. Before the
training, the data is normalized over the whole data set, so that the data is
in the value range between -1 and 1. Finally, we take advantage of the
locality of our problem and only use excerpts from the training data frames in
training. On the one hand, this saves memory, because the data is sometimes
very large and can cause problems with the GPU memory, on the other hand it
allows us to extract only relevant parts of the data. So in our example we can
only consider the data in places where there is a surface. Finally we have the
advantage to augment the data by overlapping the tiles we extract, so we can
get a lot of information from only a few frames.


The research agenda pursued in this paper was proposed by Hiroyuki Ozaki. We
are grateful to him for pointing out an important issue in economic dynamics.
This research started from an acute idea proposed by Susumu Kuwata. We
received comments from Alexander Zaslavsky and Monika Motta at the 12th AIMS
Conference on Dynamical Systems, Differential Equations and Applications. We
also received very good advice from Kenji Miyazaki at the 2018 Fall Meeting of
the Japanese Economic Association. Keiichi Morimoto also gave us sharp
comments at a workshop at Meisei University. Makoto Hanazono made a useful
suggestion at the 2020 meeting of the Japanese Society for Mathematical
Economics. Additionally, we received many comments and suggestions from
Nobusumi Sagara, Toru Maruyama, and Chaowen Yu through private discussions. We
would like to express our gratitude to all of them. Of course, the author is
responsible for all remaining errors. Finally, this work was supported by JSPS
KAKENHI Grant Number JP18K12749.

## References

[1] Acemoglu, D. 2009. Introduction to Modern Economic Growth. Princeton
University Press, Princeton.

[2] Alexandrov, A. D. 1939. Almost Everywhere Existence of the Second
Differential of a Convex Function and Some Properties of Convex Surfaces
Connected with It. Leningrad State Univ. Ann. Math. Ser. 6, 3-35.

[3] Bardi, M., Dolcetta, I, C. 2008. Optimal Control and Viscosity Solutions
of Hamilton-Jacobi-Bellman Equations. Birkhäuser, Boston.

[4] Barro, R. J., Sala-i-Martin, X. I. 2003. Economic Growth. MIT press,
Massachusetts.

[5] Benveniste, L. M., Scheinkman, J. A. 1979. On the Differentiability of the
Value Function in Dynamic Models of Economics. Econometrica, 727-732.

[6] Blanchard, O., Fischer, S. 1989. Lectures on Macroeconomics. MIT press,
Massachusetts.

[7] Coddington, E. A., Levinson, N. 1984. Theory of Ordinary Differential
Equations. McGraw-Hill.

[8] Crandall, M. G., Lions, P. L. 1983. The Viscosity Solutions of the
Hamilton-Jacobi Equations. Transactions of the American Mathematical Society
277, 1-42.

[9] Hartman, P. 1997. Ordinary Differential Equations, 2nd ed. Birkhäuser
Verlag AG.

[10] Hosoya, Y. 2019a. On the Euler Equation and the Transversality Condition
of Overtaking Optimal Solution in Macroeconomic Dynamics. IERCU discussion
paper No.320.

[11] Hosoya, Y. 2019b. The Magic of Capital. Proceedings of The Tenth
International Conference on Nonlinear Analysis and Convex Analysis, 91-108.

[12] Howard, R. 1998. Alexandrov’s Theorem on the Second Derivatives of Convex
Functions via Rademacher’s Theorem on the First Derivatives of Lipschitz
Functions. Lecture note from a functional analysis seminar at the University
of South Carolina.

[13] Ioffe, A. D., Tikhomirov, V. M. 1979. Theory of Extremal Problem. North
Holland, Amsterdam.

[14] Lions, P. L. 1982. Generalized Solutions of Hamilton-Jacobi-Bellman
Equations. Pitman Publishing, Massachusetts.

[15] Malliaris, A. G., Brock, W. A. 1988. Stochastic Methods in Economics and
Finance. North Holland.

[16] Rockafeller, R. T. 1996. Convex Analysis. Princeton University Press,
Princeton.

[17] Romer, D. 2011. Advanced Macroeconomics. McGraw-Hill.

[18] Stokey, N. L., Lucas, R. E. 1989. Recursive Methods in Macroeconomic
Dynamics. Harvard University Press, Massachusetts.