# Self-crumpling elastomers: bending induced by the drying stimulus of a
nanoparticle suspension

François Boulogne  Department of Mechanical and Aerospace Engineering,
Princeton University, Princeton, NJ 08544  Howard A. Stone  Department of
Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ
08544

###### Abstract

We report an experimental study of the drying-induced peeling of a bilayer,
consisting of an elastomeric disk coated with a suspension of nanoparticles.
We show that although capillary forces associated with the scale of the
droplet can not compete with the adhesion of the elastomer on a surface,
nevertheless large tensile stresses develop in the coating, which results in a
moment bending the bilayer. We attribute this stress to the nano-menisci in
the pores of the colloidal material and we propose a model that describes
successfully the early stage curvature of the bilayer. Thus, we show that the
peeling can be conveniently controlled by the particle size and the coating
thickness.

Many materials systems involve thin films or patches that cover substrate for
the purposes of protection, functionalization, or hierarchical design. In such
cases, it is important to be able to overcome the adhesion of the film to
substrate, as failure to do so can compromise an adhesion, damage material,
etc. In this paper, we highlight the use of a nanoparticle suspension, which,
upon drying, produces stress and bending that assist the process of
deadhesion.

To deform flexible slender bodies or thin films, two main routes have been
explored. To estimate residual or cracking stresses of consolidating thin
polymer films [1] or drying colloidal suspensions [2], the bending of a
cantilever has been used [3]. Also, locally stressed materials can fold a thin
film into a microstructured container [4], which has potential applications in
microelectronics or biomaterial encapsulation. Moreover, various stimuli have
been used to induce motion of slender structures with the use of solvent
responsive, thermoresponsive or photoactive materials [5, 6].

Capillary forces exerted by a liquid drop can also bend slender elastic
structures such as fibers [7, 8] or sheets [9, 10, 11, 12]. Considering a drop
of surface tension $\gamma$ sitting on a plate of thickness $h_{s}$ and
elastic modulus $E$, we can define an elasto-capillary length [7]
$\ell_{EC}=\sqrt{Eh_{s}^{3}/\gamma}$, which sets the critical length of the
sheet above which it will wrap around the drop. This approach of capillary
origami allows various self-folded polyhedra to be made with elastomers [10]
or metals [13].

However, to successfully achieve capillary origami, it is important to make
sure that either the adhesion of the elastomer with the surface is low or to
force the de-adhesion by an impact of the drop on the elastomer [14]. If
adhesion forces to a substrate are larger than capillary forces, then it is
impossible with a simple fluid to bend the slender elastic material. Such a
threshold is well known in peeling problems [15, 16].

Figure 1: Illustration of a typical experiment with an elastomeric disk of $7$
mm in diameter and $16$ nm particles (HS) at $\phi_{i}=0.15$. (a) Side view.
Pictures 1 to 3 show, respectively, the deposition of the colloidal droplet,
the absence of debonding by capillary forces and the removal of a fraction of
liquid to make a thin film. Pictures 4 to 6 illustrate the temporal evolution
and the corresponding pictures in (b) are the bottom view. At $t=240$ s,
enlargement of the curved bilayer (c) from a side view with a typical radius
of curvature $R=\zeta_{m}^{2}/(2\ell)$ and (d) from the bottom view. The scale
bar represents $0.5$ mm.

In this Letter, we investigate the bending of elastomeric sheets induced by
stresses developed during the drying of colloidal suspensions. We show that
where the usual elasto-capillary effect is unable to bend an elastic sheet
because of the adhesion with a substrate, debonding becomes possible by the
use of drying colloidal suspensions, which introduces the low pressure from
many small menisci in pores and leads to large tensile stresses.

## Experimental protocol

The elastomers are prepared from polydimethylsiloxane (PDMS, Sylgard 184, Dow
Corning) with a cross-linker/oligomer weight ratio of 1:10. The PDMS is spin-
coated on the top of a polystyrene petri dish lid and cured at $70^{\circ}$C
for $4$ hours and used within $3$ days. The resulting elastomer is then peeled
of the lid and carefully transferred onto a rubber sheet with a particular
attention to not stretch the film. The elastic modulus is $E_{s}=1.2$ MPa and
film thicknesses have been measured by a white light spectrometer (OceanOptics
USB2000+ used with a LS-1-LL tungsten halogen light source), confirmed by
optical profilometry (Leica DCM3D) and range from $[50,125]$
$\mathrm{\SIUnitSymbolMicro m}$. Punches are used to cut a disk in the film
with a diameter of $7$ mm (Biopsy punches, Miltex).

The elastomeric disk is then coated with a colloidal suspension (Fig. 1(a)).
In order to reproducibly make a flat coating, we follow several steps. The top
surface of the disk is made hydrophilic with a plasma treatment and the
elastomer is transferred to a transparent acrylic slide (McMaster-Carr, $4.5$
mm thickness). A drop of specified volume of colloidal suspension is deposited
with a micropipette, e.g. $20$ $\mathrm{\SIUnitSymbolMicro L}$ for a $7$ mm
disk diameter. This volume is chosen so that the drop easily spreads over the
entire surface of the disk. A given excess of liquid is removed with a second
micropipette, e.g. $16$ $\mathrm{\SIUnitSymbolMicro L}$.

The system is placed under an inverted microscope (Leica, $1.25\times$) and
pictures are recorded with a camera (Edmund Optics, CMOS monochrome USB). A
second camera (Nikon, Macro lens $105$ mm) records the time evolution of the
shape from a side view. Observations are made at relative humidity of $50\pm
5$% and temperature $23^{\circ}$C.

The colloidal suspensions we used are silica nanoparticles with different
particle diameters $2a$ and their properties are summarized in Table 1. These
suspensions can be diluted with deionized water at pH $9.5$, which is adjusted
by the addition of NaOH.

Name |  $2a$ (nm) | $\phi_{i}$ | pH
---|---|---|---
SM30 | 10 | 0.15 | 9.9
HS40 | 16 | 0.22 | 9.8
TM50 | 26 | 0.30 | 9.2
Kleb | 50 | 0.30 | 9.4
Lev | 92 | 0.30 | 8.9
Table 1: Properties of the colloidal suspensions of silica nanoparticles. The
diameter is denoted $2a$ and the initial volume fraction $\phi_{i}$. SM30,
HS40, TM50 are Ludox suspensions (purchased from Sigma-Aldrich), and Kleb and
Lev denote Klebosol 50R50 and Levasil30, respectively.

## Observations

After the deposition of the colloidal suspension on the top of the PDMS disk,
the liquid film dries from the edge, where the coating is thinner, to the
center (Fig. 1(b)) and exhibits radial crack patterns as is commonly seen in
drying colloidal droplets [17]. As the drying proceeds, tensile stresses
develop in the coating, resulting in a moment on the bilayer material (Fig.
3), which is peeled off the substrate with a macroscopic radius of curvature
$R$ as illustrated in Figs. 1(a) and (c). This phenomenon is ensured by good
adhesion of the colloidal material with the elastomer [18], which holds the
two materials constituting the bilayer together. During the delamination, in
front of the drying front, the contact area remaining between the elastomer
and the surface adopts a pattern similar to that of a flower’s petals (Fig.
1(d)), which can be attributed to the cracks resulting in mechanical
weaknesses in the bilayer. Thus, each coating strip can also bend the bilayer
in the $\theta$ direction (defined in Fig. 3) [19]. Furthermore, the
continuous peeling of the bilayer compresses its edge, which triggers a
wrinkling instability [20, 21] with a preferential location at the position of
cracks in the coating. Figure 2(b) shows the final state with different
colloidal particle sizes. Small particles induce a strong deformation of the
elastomeric sheet that could lead to several small contact points whereas
larger particles such as Klebosol have a weaker effect and form a single
contact point in the center. Finally, above a critical particle size, the
delamination is prevented and the bilayer remains flat.

The Péclet number comparing advective and diffusive effects of the particles
in the drying colloidal film can be written as
$\textrm{Pe}=6\pi\eta_{s}ah_{c}V_{E}/(k_{B}T)$, where $\eta_{s}$ is the
solvent viscosity, $V_{E}$ the evaporation rate and $k_{B}T$ the thermal
energy [22]. The drying rate of our colloidal suspension is independent of the
particle size and has a typical measured value [23] of $V_{E}\sim 10^{-8}$
m/s. This estimate is confirmed by a scaling law $V_{E}\sim h_{c}/\tau\sim
10^{-8}$ m/s where $h_{c}\sim 10$ $\mathrm{\SIUnitSymbolMicro m}$ and
$\tau\sim 100$ s is the drying timescale (Fig. 1(a)). Considering this drying
rate, we have $\textrm{Pe}\sim 10^{-2}$, which means that we can assume that
the particle distribution is homogeneous over the layer thickness.

Next, we focus on the peeling of the bilayer and we study its efficiency
through the resulting early-stage curvature for different coatings and
elastomer thicknesses. In Fig. 2(a), we present measurements of the
macroscopic curvature $R^{-1}$ of the bilayer as presented in Fig. 1(a).
Measurements are performed for a length of the debonding front
$\ell=[0.3,0.8]$ mm and we varied the particle sizes as well as the elastomer
thickness. To rationalize these results, we propose in the next paragraph a
model relating the curvature of the bilayer and the drying stress.

Figure 2: (a) Curvature $R^{-1}$ versus the particle diameter $2a$ of
colloidal suspensions at an initial volume fraction $\phi_{i}=0.15$ for
different elastomer thicknesses $h_{s}$. The disk diameter is $7$ mm. (b)
Final states (disk thickness: $50$ $\mathrm{\SIUnitSymbolMicro m}$) for
different colloidal suspensions at $\phi_{i}=0.15$.

## Model

In a classical work, Stoney studied the stress induced by a coating on a
deformable metal strip [24], assuming that both layers have the same elastic
modulus and that the coating thickness is much less than the substrate
thickness. Various extensions of this work have been reported including the
stress relief in the coating by bending [25] and the shift of the position of
the neutral axis due to different elastic moduli of the layers [26]. In our
system, the elastomer is thicker and has a lower elastic modulus than the
coating. Thus, we derive the equations under these assumptions [27].

Figure 3: Sketch of the bilayer. (a) The drying of the colloidal suspension
coated on the top of an elastomer induces tensile stresses. (b) The drying
results in a moment $M$ peeling the bilayer which remains cohesive due to a
good adhesion between colloids and the elastomer.

We assume that all quantities are independent of $\theta$ and consider the
analogous problem of a one-dimensional model for a strip. The deflection of
the material $\zeta(x)$ can be written via a Taylor series as:

$\zeta(x)=\zeta(x=0)+\partial_{x}\zeta(x=0)x+1/2\partial_{x}^{2}\zeta(x=0)x^{2}+...$

The boundary conditions are that the material is clamped at $x=0$, which leads
to $\zeta(x=0)=0$ and $\partial_{x}\zeta(x=0)=0$. Thus, we can deduce that the
beam can be considered as a portion of a circle with a radius of curvature
$R$:

$\zeta(x)=\frac{1}{2}\frac{\partial^{2}\zeta(x=0)}{\partial
x^{2}}x^{2}=\frac{x^{2}}{2R}.$ (1)

Denoting $\sigma$ the stress distribution in the beam, and $z$ the distance
across th thin film, equilibrium of the sum of axial forces and bending
moments are, respectively, expressed by:

$\displaystyle\int\sigma\,\textrm{d}S=0$ (2) $\displaystyle\int\sigma
z\,\textrm{d}S=0.$ (3)

Also, we denote the elementary cross-sectional area
$\textrm{d}S=b\,\textrm{d}z$, where $b$ is a transverse dimension, $z_{n}$ the
position of the neutral axis, $\sigma_{c}$ the stress developed in the coating
by the drying, and $h_{s}$ and $h_{c}$ the thicknesses of the substrate and
the elastomer, respectively (see Fig. 3). From equation (2), we have

$\frac{E_{s}^{\star}b}{R} \int_{0}^{h_{s}}(z_{n}-z)\,\textrm{d}z+b
\int_{h_{s}}^{h_{s}+h_{c}}\sigma_{c}\,\textrm{d}z\\\ +\frac{E_{c}^{\star}b}{R}
\int_{h_{s}}^{h_{s}+h_{c}}(z_{n}-z)\,\textrm{d}z=0,$ (4)

where $(z_{n}-z)/R$ is the strain and $E_{i}^{\star}=E_{i}/(1-\nu_{i}^{2})$
with $\nu_{i}$ the Poisson ratio of the layer $i$. The second term accounts
for the stress induced by the drying of the coating and the first and last
terms represent the response stress of each bent layer. Assuming that
$\sigma_{c}$ is independent of $z$, the integration of equation (4) determines
the position of the neutral axis $z_{n}$:

$z_{n}-h_{s}=-\frac{\sigma_{c}Rh_{c}}{E_{c}^{\star}h_{c}+E_{s}^{\star}h_{s}}+\frac{1}{2}\frac{E_{s}^{\star}h_{s}^{2}-E_{c}^{\star}h_{c}^{2}}{E_{c}^{\star}h_{c}+E_{s}^{\star}h_{s}}.$
(5)

Similarly, the integration of equation (3) with the consideration of equation
(5) leads to:

$\sigma_{c}=-\frac{E_{s}^{\star}h_{s}}{6R}P(\Sigma,\Theta),$ (6)

where we defined the dimensionless numbers
$\Sigma=E_{c}^{\star}/E_{s}^{\star}$ and $\Theta=h_{c}/h_{s}$, and the
function

$P(\Sigma,\Theta)=\frac{\Theta^{-1}+2\Sigma(2+3\Theta+2\Theta^{2})+\Sigma^{2}\Theta^{3}}{1+\Theta}.$
(7)

.

Figure 4: Dimensionless curvature ($h_{s}/R$) versus
$12\gamma/(0.15\,E_{s}^{\star}\,aP(\Sigma,\Theta))$ for different disk
thicknesses $h_{s}$. The black solid lines represent equation (9) for
$E_{c}=0.3$ GPa, without any additional prefactor. The inset represents the
radius of curvature $R$ as a function of the initial volume fraction
$\phi_{i}$ for $2a=26$ nm particles on a $7$ mm elastomer ($h_{s}=50$
$\mathrm{\SIUnitSymbolMicro m}$). The shaded areas show the domain in which
the peeling does not happen.

In the drying experiments, the air-water interface recedes into the film. We
can estimate the stress produced by the nano-menisci between colloidal
particles [22, 28]

$\sigma_{c}\simeq-\frac{2\gamma}{R_{M}}\simeq-\frac{2\gamma}{0.15a},$ (8)

where $\gamma$ is the air-water surface tension, $R_{M}\simeq 0.15a$ the
maximum radius of the meniscus and $a$ the particle radius.

Therefore, equation (8) becomes

$\frac{2\gamma}{0.15a}\simeq\frac{E_{s}^{\star}h_{s}}{6R}P(\Sigma,\Theta),$
(9)

which means that the radius of curvature $R\propto a$. Thus, the particle size
is a convenient control variable.

This prediction is compared to experimental observations on acrylic slides for
different elastomer thicknesses $h_{s}$, coating thicknesses $h_{c}$ and
particle radii $a$. Note that the coating thickness depends on the initial
volume fraction of the colloidal suspension and also on the radial position.

To quantify the effect of the coating thickness, we prepared the elastomeric
disk in the same manner as presented in the paper. These disks are glued on
glass slides by a plasma treatment to form chemical bonds. To make the disk
surface hydrophilic, a second plasma treatment is used. However, to prevent
the glass slide being hydrophilic, the glass slide not covered by the disk is
protected with tape. This step is necessary to have similar conditions to the
main experiment and to avoid the liquid to wet easily the glass slide. Thus,
it is still possible to make a large drop as shown in step 2 of Fig. 1(a).
After complete drying, we measured the shape of the coating near the disk edge
with an optical profiler. Thus, an average value of the coating thickness over
the length of the debonding front $\ell=[0.3,0.8]$ mm is deduced with an
uncertainty of $3$ $\mathrm{\SIUnitSymbolMicro m}$. We relate this thickness
to the relevant coating thickness $h_{c}$ in the delamination problem because
Fig. 1(d) shows that the peeling occurs when the material is at an advanced
step of the drying process. The coating thickness is independent of the
particle size and for $\phi_{i}=0.15$, it is about $10$
$\mathrm{\SIUnitSymbolMicro m}$.

With this approach, in Fig. 4(a), we report measurements of the radius of
curvature averaged for different experiments as a function of $a$. We must
note that the measured radius of curvature is macroscopic. As shown in figure
1(c), some deviations are observed, that can be attributed to the
inhomogeneous coating thickness.

The value of the elastic modulus of drying silica suspension is not well
established in the literature [29] due to several experimental difficulties:
heterogeneities of the material, temporal evolution due to drying, possible
effects of the drying conditions such as the geometry or the relative
humidity. We assume that $E_{c}$ is a weak function of the particle size $a$
in the studied range. We observe that our experimental results (Fig. 4) are
well described by the theoretical prediction given by equation (9) with the
fitted parameter $E_{c}=0.3$ GPa (with $\nu_{s}=0.5$ and $\nu_{c}=0.3$) in
agreement with values available in the literature [30, 23]. The collapse of
our data onto a single curve is consistent with the weak effect of the
particle size on the elastic modulus.

Moreover, we notice the limit for which the adhesion prevents the peeling of
the bilayer (Fig. 2(b) and 4). This is encountered either for large particles,
which induce smaller drying tensile stresses (equation (8)), or by a thin
coating, as illustrated by equation (6) and the inset of Fig. 4.

## Conclusion

Figure 5: Illustrations of the technique on various surfaces with a PDMS disk
with a diameter of $7$ mm and a thickness of $50$ $\mathrm{\SIUnitSymbolMicro
m}$ coated with a HS suspension at $\phi_{i}=0.15$. (a) Dust is randomly
distributed on an acrylic slide. (b) After the cleaning of the surface, more
than $90$% of the particles have been removed. The red circle shows the
initial position of the elastomeric disk. (c) Final state on a printed paper
and (d) on varnished wood (oak). Scale bars represent $1$ mm.

For a droplet deposited on the top of a thin elastomeric sheet, we have shown
that capillary forces are not sufficient to compete with adhesion and
capillary origami does not occur. Thus, in such situation for which it is
interesting to remove the sheet from the substrate, we suggest that colloidal
nanoparticles can be used to induce large tensile stresses that are able to
peel the elastomer from the substrate.

We provide a model predicting the curvature of the bilayer, which is in a good
agreement with experimental results. In particular, the model shows the effect
of the coating thickness as well as the particle size is responsible for the
tensile stress. As a result, drying tensile stresses higher than capillary
stresses can be achieved with the use of nanoparticles.

We believe that these colloidal systems can be employed in various situations
where delicate surfaces are involved such as in applications with optical and
electronic components or in restoration of photographs, painting, wallpaper
and fragile collectibles from contamination by dust, pollen, dirt, etc. [31,
32]. For illustrative purposes, we show in Fig. 5 that the bilayer studied in
this Letter can be employed to remove silica particles (typical size: $50$
$\mathrm{\SIUnitSymbolMicro m}$) randomly distributed on an acrylic slide with
a success larger than $90$% in a single step. The advantage of the method
relies on the absence of an external intervention to remove the material,
which decreases the possibility to damage the surface. We also show that the
method works on ink-printed paper or a varnished wood surface, which represent
potential applications for ancient manuscripts or antique furnitures.

## Acknowledgements

We are very grateful to F. Giorgiutti-Dauphiné and L. Pauchard for
enlightening discussions about delamination and their applications. We also
thank J.-M. Copin and M. Persson for providing Klebosol and Levasil samples
respectively. F.B. acknowledges that the research leading to these results
received funding from the People Programme (Marie Curie Actions) of the
European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant
agreement 623541. We thank the NSF for partial support via grant CBET-1132835.

## References

  * [1] K.R. Thomas and U. Steiner.  Direct stress measurements in thin polymer films.  Soft Matter, 7:7839–7842, 2011.
  * [2] H.N. Yow, M. Goikoetxea, L. Goehring, and A.F. Routh.  Effect of film thickness and particle size on cracking stresses in drying latex films.  J. Coll. and Int. Sci., 352:542–548, 2010.
  * [3] J.A. Payne.  Stress Evolution in Solidifying Coatings.  PhD thesis, 1998.
  * [4] T.G. Leong, B.R. Benson, E.K. Call, and D.H. Gracias.  Thin film stress driven self-folding of microstructured containers.  Small, 4:1605–1609, 2008.
  * [5] D.P. Holmes, M. Roche, and T. Sinha.  Bending and twisting of soft materials by non-homogenous swelling.  Soft Matter, 7:5188–5193, 2011.
  * [6] H. Meng, H. Mohamadian, M. Stubblefield, D. Jerro, S. Ibekwe, S.-S. Pang, and G. Li.  Various shape memory effects of stimuli-responsive shape memory polymers.  Smart Materials and Structures, 22:093001, 2013.
  * [7] J. Bico, B. Roman, L. Moulin, and A. Boudaoud.  Adhesion: Elastocapillary coalescence in wet hair.  Nature, 432:690–690, 2004.
  * [8] C. Duprat, S. Protiere, A.Y. Beebe, and H.A. Stone.  Wetting of flexible fibre arrays.  Nature, 482:510–513, 2012.
  * [9] J. Huang, M. Juszkiewicz, W. De Jeu, E. Cerda, T. Emrick, N. Menon, and T. Russell.  Capillary wrinkling of floating thin polymer films.  Science, 317:650–653, 2007.
  * [10] C. Py, P. Reverdy, L. Doppler, J. Bico, B. Roman, and C.N. Baroud.  Capillary origami: spontaneous wrapping of a droplet with an elastic sheet.  Phys. Rev. Lett., 98:156103, 2007.
  * [11] S. Neukirch, A. Antkowiak, and J.-J. Marigo.  The bending of an elastic beam by a liquid drop: a variational approach.  Proceedings of the Royal Society A, 469:20130066, 2013.
  * [12] A. Fargette, S. Neukirch, and A. Antkowiak.  Elastocapillary snapping: capillarity induces snap-through instabilities in small elastic beam.  Phys. Rev. Lett., 112:137802, 2014.
  * [13] T. Leong, P. Lester, T. Koh, E. Call, and D. Gracias.  Surface tension-driven self-folding polyhedra.  Langmuir, 23:8747–8751, 2007.
  * [14] A. Antkowiak, B. Audoly, C. Josserand, S. Neukirch, and M. Rivetti.  Instant fabrication and selection of folded structures using drop impact.  Proc. of the Nat. Acad. of Sci., 108:10400–10404, 2011.
  * [15] K. Kendall.  Films and Layers: Adhesion of Coatings, pages 327–351.  Springer, 2004.
  * [16] E. Hamm, P. Reis, M. LeBlanc, B. Roman, and E. Cerda.  Tearing as a test for mechanical characterization of thin adhesive films.  Nature Mater., 7:386–390, 2008.
  * [17] F. Giorgiutti-Dauphiné and L. Pauchard.  Elapsed time for crack formation during drying.  The European Physical Journal E, 37:1–7, 2014.
  * [18] L. Léger, E. Raphaël, and H. Hervet.  Surface-anchored polymer chains: their role in adhesion and friction, pages 185–225.  Springer Berlin Heidelberg, 1999.
  * [19] A.V. Zhuk, A.G. Evans, J.W. Hutchinson, and G.M. Whitesides.  The adhesion energy between polymer thin films and self-assembled monolayers.  Journal of Materials Research, 13:3555–3564, 1998.
  * [20] E. Sharon, B. Roman, M. Marder, G.-S. Shin, and H.L. Swinney.  Buckling cascades in free sheets.  Nature, 419:579–579, 2002.
  * [21] H. Aharoni and E. Sharon.  Direct observation of the temporal and spatial dynamics during crumpling.  Nature Mater., 9:993–997, 2010.
  * [22] A.F. Routh and W.B. Russel.  Horizontal drying fronts during solvent evaporation from latex films.  AIChE J., 44:2088–2098, 1998.
  * [23] F. Boulogne, L. Pauchard, and F. Giorgiutti-Dauphiné.  Effect of a non-volatile cosolvent on crack patterns induced by desiccation of a colloidal gel.  Soft Matter, 8:8505–8510, 2012.
  * [24] G.G. Stoney.  The tension of metallic films deposited by electrolysis.  Proceedings of the Royal Society of London. Series A, 82:172–175, 1909.
  * [25] E. M. Corcoran.  Determining stresses in organic coatings using plate beam deflection.  J. of Paint Tech., 41:635, 1969.
  * [26] C.-C. Chiu.  Determination of the elastic modulus and residual stresses in ceramic coatings using a strain gage.  J. Am. Ceram. Soc., 73:1999–2005, 1990.
  * [27] M. Benabdi and A.A. Roche.  Mechanical properties of thin and thick coatings applied to various substrates. Part I. An elastic analysis of residual stresses within coating materials.  Journal of Adhesion Science and Technology, 11:281–299, 1997.
  * [28] E.R. Dufresne, E.I. Corwin, N.A. Greenblatt, J. Ashmore, D.Y. Wang, A.D. Dinsmore, J.X. Cheng, X.S. Xie, J.W. Hutchinson, and D.A. Weitz.  Flow and fracture in drying nanoparticle suspensions.  Phys. Rev. Lett., 91:224501, 2003.
  * [29] L. Goehring, W. Clegg, and A. Routh.  Plasticity and fracture in drying colloidal films.  Phys. Rev. Lett., 110:024301, 2013.
  * [30] J. Zarzycki.  Critical stress intensity factors of wet gels.  J. Non-Crystalline Solids, 100:359–363, 1988.
  * [31] E. Carretti, L. Dei, and R.G. Weiss.  Soft matter and art conservation. Rheoreversible gels and beyond.  Soft Matter, 1:17–22, 2005.
  * [32] E. Carretti, M. Bonini, L. Dei, B.H. Berrie, L.V. Angelova, P. Baglioni, and R.G. Weiss.  New frontiers in materials science for art conservation: responsive gels and beyond.  Accounts of Chemical Research, 43:751–760, 2010.


# Gluformer: Transformer-Based Personalized Glucose Forecasting with
Uncertainty Quantification

###### Abstract

Deep learning models achieve state-of-the art results in predicting blood
glucose trajectories, with a wide range of architectures being proposed.
However, the adaptation of such models in clinical practice is slow, largely
due to the lack of uncertainty quantification of provided predictions. In this
work, we propose to model the future glucose trajectory conditioned on the
past as an infinite mixture of basis distributions (i.e., Gaussian, Laplace,
etc.). This change allows us to learn the uncertainty and predict more
accurately in the cases when the trajectory has a heterogeneous or multi-modal
distribution. To estimate the parameters of the predictive distribution, we
utilize the Transformer architecture. We empirically demonstrate the
superiority of our method over existing state-of-the-art techniques both in
terms of accuracy and uncertainty on the synthetic and benchmark glucose data
sets.

Index Terms—  wearable devices, time series, calibration, probabilistic
modeling

##  1 Introduction

Prediction of blood glucose values in patients with diabetes is an active area
of research [1]. From a clinical standpoint, accurate forecasting of glucose
levels can help patients take proactive actions and prevent severe
complications such as hypoglycemia or diabetic coma [2, 3]. From the
methodological perspective, the highly non-linear and non-stationary nature of
the glucose profiles makes accurate predictions difficult. To the best of our
knowledge, most current methods produce point-wise predictions without
quantifying their uncertainty, which limits clinical applicability. A
predicted glucose increase may prompt the patient to take extra insulin, which
may lead to dangerous hypoglycemia if the prediction is wrong. Complementing
predictions with uncertainty quantification will help to minimize these risks.

Since existing approaches primarily focus on maximizing the accuracy of the
predictions (rather than quantifying uncertainty), they focus on minimizing
the mean squared error (MSE) between the predictions and the observed data [4,
5, 6, 7, 8, 9, 10, 11, 12]. MSE loss implicitly corresponds to modeling future
trajectory as a function of the past plus a homogeneous Gaussian noise. If the
assumed model is correct, this approach leads to the best predictions, and
confidence intervals can be constructed by inferring the corresponding noise
variance. However, a homogeneous Gaussian noise assumption is typically
violated in practice, making the corresponding predictions and confidence
intervals unreliable.

Here we propose a new framework for modeling future glucose trajectories using
an infinite mixture model (IMM). Such formulation allows us to capture a
variety of non-Gaussian output distributions, e.g., multi-modal, skewed. We
estimate the parameters of the IMM using the Transformer architecture. While
IMM has an infinite number of parameters, we introduce an approximation step
to fit the distribution implicitly. Specifically, we inject randomness into
the computation of the network outputs and treat each output as a realization
from one of the mixture components.

Fig. 1: Predicted densities on a test sample for synthetic data set (top) and
glucose data (bottom).

To validate the capabilities of the proposed model, we consider two data sets:
1) synthetic data set 2) publicly available continuous glucose monitor (CGM)
data set [11]. We achieve state-of-the-art results for both. Figure 1 shows a
sample of predicted curves together with the visualization of predictive
distribution, demonstrating the superiority of our method in terms of the
robustness of uncertainty estimates to deviations from the Gaussian
assumption.

##  2 Related Work

Traditional machine learning. ARIMA is a well-known model for time series data
and is commonly used with CGM data [4, 5, 7, 8]. ARIMA is a non-linear auto-
regressive approach that extends the standard AR model to handle non-
stationarity by integrating the previous trends. Another popular model for
glucose predictions is random forest [6, 9, 10]. The model splits the data
into random subsamples, fitting a single prediction tree for each. The average
of the ensemble of individual trees is used for prediction. The model is
typically applied recursively by producing a single prediction at a time, but
multi-output extensions are also available.

Deep learning. As time-series data closely resemble data from the natural
language processing field, there has been an interest in applying deep
learning for multi-output time-series forecasting. Some existing state-of-the-
art approaches for glucose prediction are PolySeqMO [11] and RNN [12].
PolySeqMo models the predicted curve in terms of a polynomial expansion and
uses an RNN to learn the coefficients of the polynomial. The RNN approach of
[12] improves PolySeqMO by moving away from the polynomial expansion and
building upon the approach taken in [13], providing personalized predictions,
proposing attention for longer inputs, and introducing robust training scheme
for the model. Recently, the Temporal Fusion Transformer (TFT) has shown
promising results on a variety of time series data sets [14]. The model is
formulated using quantile regression and hence allows for uncertainty
quantification. However, unlike our model, the TFT does not provide any
measure of the distributional fit, which prevents the users from judging the
quality of the TFT uncertainty estimates.

Uncertainty estimation. The infinite mixture model (IMM) has been previously
studied in [15, 16, 17]. These approaches rely on MCMC sampling from the
posterior distribution, with latent variables explicitly being sampled from
specified prior. In contrast, we estimate the IMM parameters implicitly by
injecting latent variables (noise) into the network. This change allows us to
translate the sampling problem into the maximum-likelihood optimization, where
the only additional cost we pay is several stochastic forward passes through
the network. Since full MCMC sampling is costly to implement and diagnose for
high-dimensional data sets, our approach yields a significant computational
advantage.

Stochastic neural networks have been investigate by [18] and [19]. [18]
explore a connection between a stochastic neural network with a dropout and
$L^{2}$ regularization and the Gaussian process. Specifically, they show that
the dropout imposes a distribution on the weights of the neural network, which
can then be optimized as a variational approximation to the Gaussian kernel
parameters. Their derivations can be seen as a special case of our model if we
take the basis distribution to be Gaussian, impose $L^{2}$ constraint on the
weights, and implement the latent variables as a dropout.

##  3 Proposed model

###  3.1 Notation and background

The blood glucose prediction can be viewed as a time series forecasting
problem. Let $\mathbf{x}_{i}=(x_{ij})_{j=1}^{t}$ be a sequence of glucose
values of length $t$ with $i=\\{1,\dots,n\\}$ being the sample index. Denote
the consecutive future glucose values as $\mathbf{y}_{i}=(y_{ij})_{i=1}^{T}$,
so that $(x_{i1},\dots,x_{it},y_{i1},\dots,y_{iT})$ forms an equally-spaced
time sequence of CGM readings. Each such sequence corresponds to a unique
patient in the data set. The goal is to estimate $\mathbf{y}$ given history
$\mathbf{x}$.

To achieve this goal, previous works (Section 2) focus on building a
predictive model $f_{\bm{\theta}}:\mathbb{R}^{t}\to\mathbb{R}^{T}$ with
parameters $\bm{\theta}$. These models usually rely on the parametric
assumption
$\mathbf{y}\sim\mathcal{N}(f_{\bm{\theta}}(\mathbf{x}),\sigma^{2}I)$, which is
equivalent to the additive Gaussian assumption on the residuals. While this
assumption is sometimes explicit [13], it is often implicit through the choice
of the loss function. When the true predictive distribution,
$p(\mathbf{y}|\mathbf{x})$, is far from the Gaussian, such models are clearly
misspecified. Subsequently, both predictions and uncertainty estimates of such
models become unreliable.

###  3.2 The Infinite Mixture Model

We first consider univariate prediction and denote the scalar response
variable as $y$. We introduce latent (hidden) variables $\mathbf{z}\sim q$,
where $q$ is some distribution that allows fast sampling. Provided we have
sampled $\mathbf{z}$ from $q$, we assume that the base distribution assumption
holds conditionally on $\mathbf{z}$:

$\displaystyle y\lvert\mathbf{x},\mathbf{z}\sim
p(y;\mathbf{s}(\mathbf{x},\mathbf{z})),$

where $\mathbf{s}$ is a vector of the sufficient statistics for a distribution
$p$, e.g., mean and variance for Gaussian. Marginalizing over $\mathbf{z}$, we
get the full predictive distribution expressed as an infinite mixture model
(IMM):

$\displaystyle p(y\lvert\mathbf{x})=\int
p(y;\mathbf{s}(\mathbf{x},\mathbf{z}))q(\mathbf{z})d\mathbf{z}$

We propose to estimate the sufficient statistics $\mathbf{s}$ with a neural
network model $\hat{\mathbf{s}}\coloneqq
f_{\bm{\theta}}(\mathbf{x}|\mathbf{z})$. To approximate the intractable
predictive distribution, we use Monte Carlo integration. That is, we sample
$k$ latent $z_{j}$ independently from $q$, and average over $k$ stochastic
passes through the network:

$\mathbf{z}_{j}\sim
q,\quad\hat{\mathbf{s}}_{j}=f_{\bm{\theta}}(\mathbf{x}\lvert\mathbf{z}_{j}),\quad
p(y\lvert\mathbf{x})\approx\frac{1}{k}\sum_{j=1}^{k}p(y;\hat{\mathbf{s}}_{j}).$

The latent distribution, $q$, can be implemented by adding a random noise or
concatenating a latent noise vector to the input vector, $\mathbf{x}$.
Alternatively, the latent distribution may be embedded in the network
structure through the Bernoulli dropout layers. The choice of the base
distribution, $p$, ultimately depends on the specifics of the problem. While
the Gaussian base distribution is a good starting point due to its theoretical
properties [15], Laplace base distribution may be preferable when modeling
highly asymmetric distributions [20].

To generalize from the univariate case to the multi-output forecasting
problem, we model the multivariate predictive distribution
$p(\mathbf{y}|\mathbf{x})$ by specifying univariate IMM models as above for
each output and combining them using an appropriate copula dependence function
(e.g., independence copula, Gaussian copula etc.). Since, by Sklar’s theorem,
any multivariate distribution can be decomposed into its marginals and a
dependence function, our specification remains expressive [21]. The copula
function parameters can then also be estimated using the neural network. In
practice, we found that the independence copula, corresponding to the
independence across time assumption, was sufficient.

###  3.3 Measuring the Quality of Predictive Distribution Fit

Likelihood. To assess the quality of estimated predictive distributions, we
compute log-likelihoods for each method on the test data, where a higher value
indicates a better fit.

For existing models (Section 2), the additive Gaussian assumption on the noise
leads to the predictive distribution
$\mathbf{y}|\mathbf{x}\sim\mathcal{N}(f_{\bm{\theta}}(\mathbf{x}),\sigma^{2}I)$.
Since such models do not explicitly estimate the variance, $\sigma^{2}$,
during fitting, we estimate it using MLE and compute the average log-
likelihood as:

$\displaystyle\hat{\sigma}^{2}_{MLE}$
$\displaystyle=\frac{1}{Tn}\sum_{i=1}^{n}(\mathbf{y}_{i}-f_{\bm{\theta}}(\mathbf{x}_{i}))^{\prime}(\mathbf{y}_{i}-f_{\bm{\theta}}(\mathbf{x}_{i})),$
$\displaystyle\frac{1}{n}\log L(\bm{\theta})$
$\displaystyle=-\frac{T}{2}-\frac{T}{2}\log(2\pi\hat{\sigma}^{2}_{MLE}).$

For our model, we approximate the average log-likelihood on the test data set
based on the MC draws. That is, we estimate the average log-likelihood as:

$\displaystyle\mathbf{z}_{ij}$ $\displaystyle\sim
q,\quad\mathbf{s}_{ij}=f_{\bm{\theta}}(\mathbf{x}_{i}\lvert\mathbf{z}_{ij}),$
$\displaystyle\frac{1}{n}\log L(\bm{\theta})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\log\frac{1}{k}\sum_{j=1}^{k}p(\mathbf{y}_{i};\hat{\mathbf{s}}_{ij}).$

Calibration. To further assess the quality of our predictive distribution on
the benchmark CGM data set, we propose to evaluate model calibration. We focus
on the quality of the estimated marginal distributions for each predicted time
point. Given $n$ input sequences, our approach produces approximate predictive
distributions based on the MC draws,
$\\{\hat{p}(\mathbf{y}_{i}|\mathbf{x}_{i})\\}_{i=1}^{n}$. Based on our copula
formulation for multivariate output, we also have direct access to the
marginal distributions,
$\\{[\hat{p}(y_{ij}|\mathbf{x}_{i})]_{j=1}^{T}\\}_{i=1}^{n}$, each of which is
modeled as an IMM. Thus, we can derive approximate marginal cumulative density
functions (CDF) as:

$\hat{F}(y_{ij}|\mathbf{x}_{i})=\frac{1}{k}\sum_{i=1}^{k}F(y_{ij};\hat{\mathbf{s}}_{ij}),$

where $F$ is the CDF of the chosen univariate base distribution.

To judge the calibration of our model (how close are the estimated confidence
intervals to nominal coverage) via a diagnostic plot, we utilize the approach
described in [22]. We select a set of $\eta_{jl}\in[0,1]$ for each
$j\in\\{1,\dots,T\\}$ and a pre-specified mesh $l\in\\{1,\dots,L\\}$, and
calculate

$\displaystyle\hat{\eta}_{jl}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}\\{\hat{F}(y_{ij}|\mathbf{x}_{i})<\eta_{jl}\\}$

Thus, we obtain a set of points $\\{(\hat{\eta}_{jl},\eta_{jl})\\}_{j}$ for
each $l\in\\{1,\dots,L\\}$, where $\hat{\eta}_{jl}$ are empirical frequencies
and $\eta_{jl}$ are expected (true) frequencies. Ideally, in the case of
perfect calibration, we have that $\hat{\eta}_{jl}=\eta_{jl}$.

###  3.4 Estimation

Transformer Architecture. The Transformer model [23] consists of a series of
attention heads and fully connected layers. To overcome the computational cost
of modeling longer prediction intervals, we use a modification of the
Transformer architecture that introduces convolutional layers in-between
attention blocks and is able to provide multi-output predictions in a single
pass [24]. We use the model to predict the parameters of each of the marginal
IMM’s, where we consider each stochastic pass through the model to be an MC
draw.

Feature Embedding. As outlined in Section 4.1, the CGM data usually only
contains glucose readings together with the timestamps. It is reasonable to
assume that each subject adheres to a daily schedule that affects his or her
blood glucose levels. Therefore, we utilize time features to enhance our
predictions. Also, we add positional and subject encoding.

Fig. 2: Calibration curves depicting $\\{(\hat{\eta}_{jl},\eta_{jl})\\}_{j}$
for $l=\\{1,\dots,12\\}$ for proposed model on the benchmark CGM data. Perfect
calibration would be represented by the $45^{\circ}$ line.

Model | IMM | Gaussian
---|---|---
APE | 28.63 | 32.68
RMSE | 0.1355 | 0.1415
Log-likelihood | -0.65 | -4.21

Table 1: APE, RMSE, and likelihood computed on the test set for synthetic
data.

Model | Full | Event | Hypo | Hyper | Likelihood
---|---|---|---|---|---
ARIMA | 9.85 / 17.65 | 8.91 / 19.86 | 19.94 / 14.53 | 8.51 / 22.17 | -14.93
RF:Rec | 9.04 / 17.15 | 8.97 / 20.36 | 18.84 / 12.43 | 8.68 / 23.41 | -14.58
RF:MO | 10.22 / 18.27 | 8.61 / 19.90 | 21.64 / 17.36 | 7.99 / 21.58 | -15.34
PolySeqMO | 8.55 / 15.68 | 8.27 / 18.81 | 22.86 / 21.87 | 6.77 / 18.30 | -15.61
RNN | 8.17 / 15.67 | 8.29 / 19.37 | 18.72 / 16.26 | 6.99 / 19.22 | -13.50
TFT | 7.80 / 15.78 | 8.03 / 18.23 | 16.23 / 14.62 | 6.87 / 18.98 | –
Our: IMM |  7.78 / 15.40 | 7.89 / 17.85 |  15.75 / 14.03 | 7.08 / 19.58 | -2.67
Our: Gaussian | 7.82 / 14.73 | 7.93 / 18.62 | 18.28 / 13.55 | 6.81 / 18.94 | -12.82

Table 2: APE/RMSE for 60-minute prediction window (Full), hypoglycemia,
hyperglycemia, and event (hypo-or hyperglycemia), and model log-likelihood on
test data.

##  4 Performance Evaluation

###  4.1 Data

Synthetic. We consider a sample of $2000$ time series for training, $100$ for
validation, and $100$ for testing. We generate the data using a mixture of
Gaussian Processes. Each series has the same range of values in the time
interval $[-4,0)$, but in the interval $[0,3]$ takes one of two trajectories:
increase or decrease. The resulting distribution has a varying number of
modes: one mode on the interval $[-4,0)$ and two modes on $[0,3]$. The
training curves are depicted in grey in Figure 1.

Benchmark CGM. We use a publicly available CGM data set [11], which contains
information on glucose levels of $38$ subjects tracked continuously throughout
multiple disjoint intervals with the measurement frequency of $5$ minutes.
Similar to [12], we remove periods of drastic fluctuations where the
subsequent CGM readings change by more than $40~{}\text{mg/dL}$. To be
compatible with the analysis in [12], we do not utilize any interpolation
techniques for the missing data. The resulting data consists of $399,302$
observations, with an average of 30 (uninterrupted) sequences per subject of
length between 200 and 400 observations. We split the data by sequences into
the train, validation, and test sets in $20:1:1$ proportion.

###  4.2 Model Implementation

We implement our method in PyTorch [25], and use a single NVIDIA RTX 2080 Ti
GPU for all experiments. For ablation purposes, we train the model using (i)
proposed loss corresponding to the IMM; (ii) mean squared error (MSE)
corresponding to the Gaussian assumption. For IMM, we implement the latent
variables by introducing the dropout and assuming independence across time.
The full model implementation is available at github.com/gluformer. For the
synthetic data, we set the encoder length to 4 and the prediction to 2. For
the glucose data, we set the encoder length to 180 (36 hours) and the
prediction length to 12 (1 hour).

###  4.3 Results

We measure prediction accuracy using root mean squared error (RMSE) and
average percentage error (APE) across prediction intervals. We report the
final APE and RMSE as medians of the individual prediction errors across
samples. We measure the quality of fit as in Section 3.3.

Synthetic. We compare the proposed IMM to the same architecture trained with
the Gaussian assumption. We report our results in Table 1 and plot a sample of
predictions in Figure 1, displaying the predicted densities for the two
models. The proposed IMM is more accurate and correctly captures the multi-
modality in the predicted densities across time points, providing a better
probabilistic fit to the data.

Benchmark CGM. We report both the error on the full data set as well as errors
for the periods of hypo- and hyperglycemia [26]. We define a period of
hypoglycemia if the glucose level drops below $70~{}\text{mg/dL}$ in the
forecast window and a period of hyperglycemia if the glucose level rises above
$180~{}\text{mg/dL}$. We also measure the model’s accuracy during the
incidence of either of the events. Table 2 summarize the results on the
benchmark CGM data for the 60 minutes prediction window. The proposed IMM
outperforms the others 3 out of 4 times in APE and 2 out of 4 times in RMSE.
In particular, IMM achieves state-of-the-art $7.78$ APE for the most difficult
$60$ minute forecast. Table 2 also reports average likelihoods on test data
for all methods, except TFT for which it is not available. Our model achieves
$6\times$ better average log-likelihood than the previous methods. We believe
this to be indicative of the fact that the true predictive distribution of
$\mathbf{y}|\mathbf{x}$ is not similar to Gaussian (likely multi-modal). We
report our calibration plots for each of the $12$ prediction points in the
future (the longest prediction length) in Figure 2. When plotted, a perfectly
calibrated model should correspond to a straight line. We find that the model
observed confidence levels are close to the expected levels for all the
points, highlighting the model’s excellent performance.

##  5 Conclusion

We propose a novel approach to quantify uncertainty in multi-output time
series predictions based on the infinite mixture model. We utilize an
efficient Transformer architecture for estimation. When applied to public CGM
data, our model outperforms existing approaches (ARIMA random forest, RNN, and
Transformer) in terms of the achieved log-likelihood on the test data and
provides competitive accuracy for various forecast horizons and dangerous
event scenarios.

## References

  * [1] Silvia Oviedo, Josep Vehí, Remei Calm, and Joaquim Armengol,  “A review of personalized blood glucose prediction strategies for t1dm patients,”  International journal for numerical methods in biomedical engineering, vol. 33, no. 6, pp. e2833, 2017.
  * [2] Hong Sun, Pouya Saeedi, Suvi Karuranga, Moritz Pinkepank, Katherine Ogurtsova, Bruce B Duncan, Caroline Stein, Abdul Basit, Juliana CN Chan, Jean Claude Mbanya, et al.,  “Idf diabetes atlas: Global, regional and country-level diabetes prevalence estimates for 2021 and projections for 2045,”  Diabetes research and clinical practice, p. 109119, 2021.
  * [3] Steven Broll, Jacek Urbanek, David Buchanan, Elizabeth Chun, John Muschelli, Naresh M Punjabi, and Irina Gaynanova,  “Interpreting blood glucose data with r package iglu,”  PloS one, vol. 16, no. 4, pp. e0248560, 2021.
  * [4] Meriyan Eren-Oruklu, Ali Cinar, and Lauretta Quinn,  “Hypoglycemia prediction with subject-specific recursive time-series models,” 2010.
  * [5] Dimitri Boiroux, Anne Katrine Duun-Henriksen, Signe Schmidt, Kirsten Nørgaard, Sten Madsbad, Ole Skyggebjerg, Peter Ruhdal Jensen, Niels Kjølstad Poulsen, Henrik Madsen, and John Bagterp Jørgensen,  “Overnight control of blood glucose in people with type 1 diabetes,”  IFAC Proceedings Volumes, vol. 45, no. 18, pp. 73–78, 2012.
  * [6] Bharath Sudharsan, Malinda Peeples, and Mansur Shomali,  “Hypoglycemia prediction using machine learning models for patients with type 2 diabetes,”  Journal of diabetes science and technology, vol. 9, no. 1, pp. 86–90, 2014.
  * [7] Mwaffaq Otoom, Hussam Alshraideh, Hisham M Almasaeid, Diego López-de Ipiña, and José Bravo,  “Real-time statistical modeling of blood sugar,”  Journal of medical systems, vol. 39, no. 10, pp. 1–6, 2015.
  * [8] Carlo Novara, N Mohammad Pour, Tyrone Vincent, and Giorgio Grassi,  “A nonlinear blind identification approach to modeling of diabetic patients,”  IEEE Transactions on Control Systems Technology, vol. 24, no. 3, pp. 1092–1100, 2015.
  * [9] Eleni I Georga, Vasilios C Protopappas, Demosthenes Polyzos, and Dimitrios I Fotiadis,  “Evaluation of short-term predictors of glucose concentration in type 1 diabetes combining feature ranking with regression models,”  Medical & biological engineering & computing, vol. 53, no. 12, pp. 1305–1318, 2015.
  * [10] J Ignacio Hidalgo, J Manuel Colmenar, Gabriel Kronberger, Stephan M Winkler, Oscar Garnica, and Juan Lanchares,  “Data based prediction of blood glucose concentrations using evolutionary methods,”  Journal of medical systems, vol. 41, no. 9, pp. 1–20, 2017.
  * [11] Ian Fox, Lynn Ang, Mamta Jaiswal, Rodica Pop-Busui, and Jenna Wiens,  “Deep Multi-Output Forecasting: Learning to Accurately Predict Blood Glucose Trajectories,”  ACM SIGKDD, pp. 1387–1395, July 2018.
  * [12] Mohammadreza Armandpour, Brian Kidd, Yu Du, and Jianhua Z. Huang,  “Deep Personalized Glucose Level Forecasting Using Attention-based Recurrent Neural Networks,”  IJCNN, June 2021.
  * [13] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski,  “Deepar: Probabilistic forecasting with autoregressive recurrent networks,”  International Journal of Forecasting, vol. 36, no. 3, pp. 1181–1191, 2020.
  * [14] Bryan Lim, Sercan O Arik, Nicolas Loeff, and Tomas Pfister,  “Temporal fusion transformers for interpretable multi-horizon time series forecasting,”  International Journal of Forecasting, 2019.
  * [15] Carl Rasmussen,  “The infinite gaussian mixture model,”  NeurIPS, vol. 12, 1999.
  * [16] Ian Porteous, Alexander T Ihler, Padhraic Smyth, and Max Welling,  “Gibbs sampling for (coupled) infinite mixture models in the stick breaking representation,”  arXiv preprint arXiv:1206.6845, 2012.
  * [17] Halid Z Yerebakan, Bartek Rajwa, and Murat Dundar,  “The infinite mixture of infinite gaussian mixtures,”  NeurIPS, vol. 27, 2014.
  * [18] Yarin Gal and Zoubin Ghahramani,  “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,”  in ICML, Maria Florina Balcan and Kilian Q. Weinberger, Eds., New York, New York, USA, 20–22 Jun 2016, vol. 48 of PMLR, pp. 1050–1059, PMLR.
  * [19] Stefan Depeweg, José Miguel Hernández-Lobato, Finale Doshi-Velez, and Steffen Udluft,  “Learning and policy search in stochastic dynamical systems with bayesian neural networks,”  in ICLR, 2017.
  * [20] Axel Brando, Jose A Rodriguez, Jordi Vitria, and Alberto Rubio Muñoz,  “Modelling heterogeneous distributions with an uncountable mixture of asymmetric laplacians,”  NeurIPS, vol. 32, 2019.
  * [21] M Sklar,  “Fonctions de repartition an dimensions et leurs marges,”  Publ. inst. statist. univ. Paris, vol. 8, pp. 229–231, 1959.
  * [22] Volodymyr Kuleshov, Nathan Fenner, and Stefano Ermon,  “Accurate uncertainties for deep learning using calibrated regression,”  in ICML. PMLR, 2018, pp. 2796–2804.
  * [23] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin,  “Attention is all you need,”  in NeurIPS, 2017, pp. 5998–6008.
  * [24] Haoyi Zhou, Shanghang Zhang, Jieqi Peng, Shuai Zhang, Jianxin Li, Hui Xiong, and Wancai Zhang,  “Informer: Beyond Efficient Transformer for Long Sequence Time-Series Forecasting,”  arXiv:2012.07436 [cs], Mar. 2021.
  * [25] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.,  “Pytorch: An imperative style, high-performance deep learning library,”  NeurIPS, vol. 32, 2019.
  * [26] Kamuran Turksoy, Elif S Bayrak, Lauretta Quinn, Elizabeth Littlejohn, Derrick Rollins, and Ali Cinar,  “Hypoglycemia early alarm systems based on multivariable models,”  Industrial & engineering chemistry research, vol. 52, no. 35, pp. 12329–12336, 2013.