# Nonlinear and Nonlocal Elasticity in Coarse-Grained Differential-Tension
Models of Epithelia

Pierre A. Haas<EMAIL_ADDRESS>Raymond E. Goldstein
<EMAIL_ADDRESS>Department of Applied Mathematics and
Theoretical Physics, Centre for Mathematical Sciences,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom

###### Abstract

The shapes of epithelial tissues result from a complex interplay of
contractile forces in the cytoskeleta of the cells in the tissue, and adhesion
forces between them. A host of discrete, cell-based models describe these
forces by assigning different surface tensions to the apical, basal, and
lateral sides of the cells. These differential-tension models have been used
to describe the deformations of epithelia in different living systems, but the
underlying continuum mechanics at the scale of the epithelium are still
unclear. Here, we derive a continuum theory for a simple differential-tension
model of a two-dimensional epithelium and study the buckling of this
epithelium under imposed compression. The analysis reveals how the cell-level
properties encoded in the differential-tension model lead to linear, nonlinear
as well as nonlocal elastic behavior at the continuum level.

##  I Introduction

Intercellular adhesion proteins and cortical actin networks are well
established as regulators of cell surface mechanics, and hence of the
deformations of epithelia during morphogenesis Lecuit and Lenne (2007). Ever
since the seminal work of Odell et al. Odell  _et al._ (1981), these cellular
components have therefore underlain mathematical models of epithelia Rauzi
_et al._ (2013); Fletcher  _et al._ (2016). One large class of such models are
differential-tension models Derganc  _et al._ (2009); Brezavšček  _et al._
(2012); Krajnc  _et al._ (2013); Hannezo  _et al._ (2014); Rauzi  _et al._
(2015); Krajnc and Ziherl (2015); Štorgel  _et al._ (2016); Krajnc  _et al._
(2018), in which cell polarity, cell-cell adhesion properties and the
actomyosin network induce different surface tensions in different sides of the
discrete cells. Coupling the mechanical models describing epithelial
deformations to models of the intracellular biochemistry is a key challenge in
the field Howard  _et al._ (2011), but some progress has recently been made by
coupling the differential-tension model of Ref. Hannezo  _et al._ (2014) to
the diffusion of a ‘mechanogen’ that induces contractility Dasbiswas  _et al._
(2018).

While the differential-tension models can quantitatively reproduce the
morphology of epithelial folds in many different living systems Štorgel  _et
al._ (2016), it is likely that, in general multilayered epithelia, the
formation of epithelial folds must be ascribed to a combination of these
intra-epithelial stresses and differential growth of different parts of the
tissue. Models based on the latter only have for example been invoked to
describe, at the scale of the epithelium, the formation of cortical
convolutions in the brain Richman  _et al._ (1975); Tallinen  _et al._ (2013);
Manyuhina  _et al._ (2014); Budday  _et al._ (2015); Tallinen  _et al._
(2016); Lejeune  _et al._ (2016) and of the intestinal villi Hannezo  _et al._
(2011); Savin  _et al._ (2011); Shyer  _et al._ (2013), the ‘fication’ of
which lends itself to the pun that gave Ref. Shyer  _et al._ (2013) its title.
The coarse-grained limit of the differential-tension models at this scale is
less well-studied, however, and this continuum limit is the topic of this
paper.

We shall focus on the most basic setup of these differential-tension models
Derganc  _et al._ (2009); Krajnc  _et al._ (2013); Krajnc and Ziherl (2015) in
an epithelial monolayer: the apical, basal, and lateral sides of the cells
have respective areas $A_{\mathrm{a}}$, $A_{\mathrm{b}}$, and $A_{\ell}$. The
internal state of the cells induces different surface tensions
$\Gamma_{\mathrm{a}}$, $\Gamma_{\mathrm{b}}$, and $\Gamma_{\ell}$ in the
apical, basal, and lateral sides of the cells, respectively. The energy of a
single cell therefore reads

$\displaystyle
E=\Gamma_{\mathrm{a}}A_{\mathrm{a}}+\Gamma_{\mathrm{b}}A_{\mathrm{b}}+\dfrac{\Gamma_{\ell}}{2}A_{\ell},$
(1)

where the factor of $1/2$ has been introduced for mere convenience Krajnc  _et
al._ (2013); Krajnc and Ziherl (2015). The theory can be extended to
incorporate additional physics such as a basement membrane Štorgel  _et al._
(2016) or a confining vitelline membrane Brezavšček  _et al._ (2012). In a
full three-dimensional setup, this leads to the study of the shapes of prism-
shaped cells Hannezo  _et al._ (2014); Krajnc  _et al._ (2018). Here, we shall
restrict to the two-dimensional setup Krajnc  _et al._ (2013); Krajnc and
Ziherl (2015) of an epithelium consisting of isosceles trapezoidal cells of
parallel apical and basal sides of lengths $L_{\mathrm{a}}$ and
$L_{\mathrm{b}}$. These trapezoids are joined up across their lateral sides,
which have equal length $L_{\mathrm{\ell}}$ (Fig. 1). Since there are two such
lateral sides, the cell energy (1) reduces to

$\displaystyle
E=\Gamma_{\mathrm{a}}L_{\mathrm{a}}+\Gamma_{\mathrm{b}}L_{\mathrm{b}}+\Gamma_{\ell}L_{\ell},$
(2)

per unit extent in the third dimension Krajnc  _et al._ (2013); Krajnc and
Ziherl (2015). A different two-dimensional limit is obtained by averaging over
the thickness of the cell sheet and describing in-plane deformations only.
Such models, termed area- and perimeter-elasticity models, have been studied
extensively Hufnagel  _et al._ (2007); Farhadifar  _et al._ (2007);
Hilgenfeldt  _et al._ (2008); Bi  _et al._ (2015, 2016).

Figure 1: Model epithelium. Isosceles trapezoidal cells of apical and basal
bases $L_{\mathrm{a}}$ and $L_{\mathrm{b}}$ are connected along their lateral
sides, which have length $L_{\ell}$. The continuous line and shaded area
provide a cartoon of the continuum limit.

The simplest problem in the mechanics of elastic rods is their Euler buckling
under applied forces Landau and Lifshitz (1970); it is therefore meet to ask
how the buckling behavior of an active material such as this model epithelium
differs from that of an elastic rod. This problem was considered in Ref.
Krajnc and Ziherl (2015), where the continuum limit of Eq. (2) was mapped to
Euler’s Elastica equation *[][App.A; pp.546--557;Chap.6.4; pp.184--
192andApp.D; pp.571--581.]audolypomeau. These calculations, complementing
simulations of the discrete model in Ref. Krajnc  _et al._ (2013), were
phrased in terms of spontaneous buckling of the epithelium, but an additional
compressive or extensile force is required to produce these deformations,
making the analysis of Ref. Krajnc and Ziherl (2015) more appropriate to the
present setup of buckling under imposed forces. Moreover, the analysis of Ref.
Krajnc and Ziherl (2015) is not completely consistent with the discrete model,
since it does not impose the condition that the trapezoidal cells match up
exactly along their lateral sides.

Figure 2: Cell Geometry. (a) Geometry of a single isosceles trapezoidal cell
of mean base $K$ and height $L$, and sidelengths
$L_{\mathrm{a}},L_{\mathrm{b}},L_{\ell}$, the lateral sides being at an angle
$2\phi$ to each other. (b) Definition of the tangent angle $\psi$ of the
midline below the horizontal. The geometry of contiguous cells defines the
relation between $\phi$ and $\psi$, as expressed in Eq. (11).

In this paper, we perform a consistent asymptotic expansion of the discrete
geometry of this model, revealing nonlinear and nonlocal elastic terms in the
continuum limit. We then analyze the buckling behavior of the continuum model
under imposed compression analytically and numerically.

##  II Continuum Model

###  II.1 Single-cell Energy

To obtain the energy of a single cell, we express the sidelengths of the
trapezoidal cells in terms of their mean base $K$, their height $L$, and the
angle $2\phi$ that their lateral sides make with each other (Fig. 2a):

$\displaystyle L_{\mathrm{a}}$ $\displaystyle=K+L_{\ell}\sin{\phi},$ (3a)
$\displaystyle L_{\mathrm{b}}$ $\displaystyle=K-L_{\ell}\sin{\phi},$ (3b)
$\displaystyle L$ $\displaystyle=L_{\ell}\cos{\phi}.$ (3c) Incompressibility
implies the cell area conservation constraint $\displaystyle
A_{\mathrm{c}}=KL.$ (3d)

Upon eliminating $K$ and $L$ using these relations, the energy of a single
cell is expressed, from Eq. (2) and as a function of $L_{\ell}$ and $\phi$, as

$\displaystyle
E=(\Gamma_{\mathrm{a}}+\Gamma_{\mathrm{b}})\dfrac{A_{\mathrm{c}}}{L_{\ell}}\sec{\phi}+(\Gamma_{\mathrm{a}}-\Gamma_{\mathrm{b}})L_{\ell}\sin{\phi}+\Gamma_{\ell}L_{\ell}.$
(4)

####  II.1.1 Non-dimensionalization

We non-dimensionalize this expression by scaling lengths with the square root
of the cell area and thus define the non-dimensional length of the lateral
sides of the trapezoidal cells,
$\smash{\lambda=L_{\ell}/A_{\mathrm{c}}^{1/2}}$. We further set
$\smash{\ell=L/A_{\mathrm{c}}^{1/2}}$ and
$\smash{\kappa=K/A_{\mathrm{c}}^{1/2}}$. Finally, following Ref. Krajnc and
Ziherl (2015), we introduce the parameters

$\displaystyle\alpha=\dfrac{\Gamma_{\mathrm{a}}}{\Gamma_{\ell}},$
$\displaystyle\beta=\dfrac{\Gamma_{\mathrm{b}}}{\Gamma_{\ell}},$
$\displaystyle\ell_{0}=\sqrt{\alpha+\beta},$
$\displaystyle\delta=\alpha-\beta.$ (5)

We note that $\ell_{0}$ is the (uniform and non-dimensionalized) thickness of
the epithelium in the flat configuration. Hence $s_{0}=1/\ell_{0}$ is the
(non-dimensionalized) width of a single cell (i.e. its arclength in the flat
configuration). The non-dimensional energy
$\smash{e=E/\Gamma_{\ell}A_{\mathrm{c}}^{1/2}}$ of a single cell is therefore

$\displaystyle
e=\dfrac{\ell_{0}^{2}}{\lambda}\sec{\phi}+\lambda\bigl{(}\delta\sin{\phi}+1\bigr{)}.$
(6)

Without loss of generality, we assume that $\delta>0$ in what follows, so, for
a single cell, $\phi<0$ is energetically favorable.

####  II.1.2 Transition to Constricted Cells

The transition to constricted triangular cells is a geometric singularity in
the discrete model. These triangular cells arize as limiting cases of the
trapezoidal cells when $L_{\text{a}}=0$ or $L_{\text{b}}=0$. Using Eq. (3),
the conditions $L_{\text{a}},L_{\text{b}}\geqslant 0$ reduce to

$\displaystyle\lambda^{2}\leqslant\dfrac{1}{\cos{\phi}\left|\sin{\phi}\right|}=\dfrac{2}{\left|\sin{2\phi}\right|}.$
(7)

###  II.2 Energy of an Epithelium

In the continuum limit, we take $\phi$ to be a function of the arclength $s$
of the midline of the undeformed, flat epithelium. Summing over all the cells,
we obtain the non-dimensional energy $\mathcal{E}$ of the epithelium,

$\displaystyle\mathcal{E}=\int{e(\phi)\,\ell_{0}\,\mathrm{d}s}=\ell_{0}^{2}\int{\left(\dfrac{\ell_{0}}{\lambda}\sec{\phi}+\dfrac{\lambda}{\ell_{0}}\bigl{(}1+\delta\sin{\phi}\bigr{)}\right)\,\mathrm{d}s},$
(8)

By imposing an energy density equal to $1/s_{0}=\ell_{0}$ times the (non-
dimensional) energy of a single cell and integrating with respect to the
reference arclength in this manner, we have imposed local cell area
conservation 111In this respect, the present analysis differs from that of
Ref. Krajnc and Ziherl (2015), where the energy was expressed as an integral
with respect to the arclength $S$ in the deformed configuration, related in
fact to $s$ by $\ell\,\mathrm{d}S=\ell_{0}\,\mathrm{d}s$. In Ref. Krajnc and
Ziherl (2015), a weaker global area conservation constraint for the cell sheet
had therefore to be imposed separately..

The boundary conditions are most naturally expressed in terms of the angle
$\psi$ of the deformed midline of the epithelium _below_ the horizontal (Fig.
2b). We therefore express the energy in terms of $\psi$. This is usefully done
in the scaling limit $\ell_{0}\gg 1$ of a columnar epithelium, as explained
next.

####  II.2.1 Asymptotic Expansion

We make two further scaling assumptions:

$\displaystyle\lambda=\mathcal{O}(\ell_{0}),$
$\displaystyle\phi=\mathcal{O}\bigl{(}\ell_{0}^{-2}\bigr{)}.$ (9)

These scalings correspond to the regime $\ell\sim\ell_{0}$ and $\phi\ll 1$,
where the cells deform but slightly from their equilibrium configuration. In
this limit, Eq. (3c) implies that $\lambda\sim\ell\sim\ell_{0}$. Further, area
conservation (3d) requires that $\kappa\sim 1/\ell_{0}$, and thus, from Eqs.
(3a,3b), we must have $\phi\lesssim\kappa/\lambda\sim 1/\ell_{0}^{2}$. The
second scaling thus corresponds to the largest deformations allowed. We
therefore introduce the parameter

$\displaystyle\Lambda=\lambda/\ell_{0}=\mathcal{O}(1).$ (10)

We are now set up to relate $\phi$ and $\psi$, for which purpose we use the
geometric relation

$\displaystyle\psi(s+ks_{0})-\psi(s)$
$\displaystyle=\phi(s)+2\phi(s+s_{0})+\cdots$ $\displaystyle\hskip
14.22636pt+2\phi\bigl{(}s\\!+\\!(k\\!-\\!1)s_{0}\bigr{)}+\phi(s\\!+\\!ks_{0}),$
(11)

valid for any positive integer $k$, as sketched in Fig. 2b. In Appendix A, we
show that, with our scaling assumptions, the continuum limit of this relation
is

$\displaystyle\psi^{\prime}(s)=\sum_{m=0}^{\infty}{\dfrac{2\mathcal{B}_{2m}}{(2m)!}\dfrac{\phi^{(2m)}(s)}{\ell_{0}^{2m-1}}},$
(12)

wherein $\mathcal{B}_{0}=1,\mathcal{B}_{1}=-\frac{1}{2},\dots$ are the
Bernoulli numbers (of the first kind) Spanier and Oldham (1987). The next step
is to invert this series, to express $\phi$ in terms of the derivatives of
$\psi$. While we are not aware of any explicit expression for the coefficients
of the inverted series, it is straightforward to invert the series order-by-
order by substituting back and forth, and thus obtain

$\displaystyle\phi(s)=\dfrac{\psi^{\prime}(s)}{2\ell_{0}}-\dfrac{\psi^{\prime\prime\prime}(s)}{24\ell_{0}^{3}}+\dfrac{\psi^{(\text{v})}(s)}{240\ell_{0}^{5}}+\cdots,$
(13)

where dashes denote differentiation with respect to $s$. In Ref. Krajnc and
Ziherl (2015), only the first term of this expansion was obtained. Inclusion
of the second term will enable us to analyze the buckling behavior of the
epithelium in what follows.

####  II.2.2 Shape Equations for the Buckled Epithelium

We describe the shape of the buckled epithelium by the coordinates
$\bigl{(}x(s),y(s)\bigr{)}$ of the centreline of the epithelium, defined by
the axes in Fig. 2b. To derive the continuum equations describing the
centreline, we begin by projecting the discrete geometry onto the axes,

$\displaystyle x(s\\!+\\!s_{0})-x(s)$
$\displaystyle=\dfrac{1}{2}\Bigl{(}\kappa(s)\cos{\psi(s)}+\kappa(s\\!+\\!s_{0})\cos{\psi(s\\!+\\!s_{0})}\Bigr{)},$
(14a) $\displaystyle y(s\\!+\\!s_{0})-y(s)$
$\displaystyle=\dfrac{1}{2}\Bigl{(}\kappa(s)\sin{\psi(s)}+\kappa(s\\!+\\!s_{0})\sin{\psi(s\\!+\\!s_{0})}\Bigr{)}.$
(14b)

Using $\kappa(s)=(\ell_{0}\Lambda)^{-1}\sec{\phi(s)}$ and expanding these
equations order-by-order in inverse powers of $\ell_{0}$ using Eq. (13), we
obtain, after a considerable amount of algebra 222These and other expansions
in this paper were carried out with the help of Mathematica (Wolfram, Inc.) to
facilitate the manipulation of complicated algebraic expressions.,

$\displaystyle\Lambda\dfrac{\mathrm{d}x}{\mathrm{d}s}=f\cos{\psi}-g\sin{\psi},$
$\displaystyle\Lambda\dfrac{\mathrm{d}y}{\mathrm{d}s}=f\sin{\psi}+g\cos{\psi},$
(15)

where

$\displaystyle f$
$\displaystyle=1+\dfrac{{\psi^{\prime}}^{2}}{24\ell_{0}^{2}}+\dfrac{7{\psi^{\prime}}^{4}+144{\psi^{\prime\prime}}^{2}+32\psi^{\prime}\psi^{\prime\prime\prime}}{5760\ell_{0}^{4}}+\mathcal{O}\bigl{(}\ell_{0}^{-6}\bigr{)},$
(16a) $\displaystyle g$
$\displaystyle=\dfrac{\psi^{\prime\prime}}{12\ell_{0}^{2}}+\dfrac{87{\psi^{\prime}}^{2}\psi^{\prime\prime}-2\psi^{\text{(iv)}}}{1440\ell_{0}^{4}}+\mathcal{O}\bigl{(}\ell_{0}^{-6}\bigr{)}.$
(16b)

Integrating these differential equations yields the shape of the buckled
epithelium. Deviations from the ‘standard’ values $f=1$, $g=0$ arise at order
$\smash{\mathcal{O}\bigl{(}\ell_{0}^{-2}\bigr{)}}$.

####  II.2.3 Derivation of the Governing Equation

We shall seek to describe buckled configurations of an epithelium of
undeformed length $2\Sigma\gg s_{0}$. We change variables by introducing
$\sigma=s/\Sigma$, use dots to denote differentiation with respect to
$\sigma$, and define

$\displaystyle\Xi=\ell_{0}\Sigma\gtrsim\mathcal{O}(1).$ (17)

Since $\Xi=\Sigma/s_{0}$, the number of cells in the epithelium is simply
$N=2\Xi$.

We shall seek buckled solutions with clamped boundary conditions and a
prescribed relative compression $D$, so that

$\displaystyle x(2)-x(0)=2(1-D),$ $\displaystyle y(2)=y(0),$ (18)

where the coordinates are now expressed relative to the scaled arclength
$\sigma$. We shall restrict to symmetrically buckled configurations for which
$\psi(\sigma)=-\psi(2-\sigma)$. The second condition above is then satisfied.
We may further reduce the solution to the range $0\leqslant\sigma\leqslant 1$,
with the condition of prescribed compression reading $x(1)-x(0)=1-D$. To
minimize the energy of the epithelium at this imposed displacement, we
therefore consider the Lagrangian

$\displaystyle\mathcal{L}$
$\displaystyle=\int_{0}^{1}{\left(\dfrac{\sec{\phi(\sigma)}}{\Lambda}+\Lambda\bigl{(}1+\delta\sin{\phi(\sigma)}\bigr{)}\right)\mathrm{d}\sigma}$
$\displaystyle\hskip
128.0374pt+\dfrac{\mu}{\Sigma}\int_{0}^{1}{\dot{x}(\sigma)\,\mathrm{d}\sigma},$
(19)

where the Lagrange multiplier $\mu$ imposes the displacement condition and has
the interpretation of a horizontal, compressive force.

Upon substituting for $\phi$ using Eq. (13), expanding in inverse powers of
$\Xi$, discarding terms that vanish upon integration, and integrating by
parts, we find

$\displaystyle\mathcal{L}=\int_{0}^{1}{\left[\Lambda+\dfrac{1}{\Lambda}+\dfrac{\dot{\psi}^{2}}{8\Lambda\Xi^{2}}-\dfrac{\delta\Lambda\dot{\psi}^{3}}{48\Xi^{3}}+\dfrac{5\dot{\psi}^{4}}{384\Lambda\Xi^{4}}+\dfrac{\ddot{\psi}^{2}}{48\Lambda\Xi^{4}}+\dfrac{\mu\cos{\psi}}{\Lambda}\left(1+\dfrac{\dot{\psi}^{2}}{8\Xi^{2}}+\dfrac{41\dot{\psi}^{4}}{1920\Xi^{4}}+\dfrac{\ddot{\psi}^{2}}{40\Xi^{4}}+\dfrac{\dot{\psi}\dddot{\psi}}{240\Xi^{4}}\right)+\mathcal{O}\bigl{(}\Xi^{-5}\bigr{)}\right]\mathrm{d}\sigma}.$
(20)

To analyze the dependence of the energy on the differential tension $\delta$,
we must go beyond lowest order 333A non-trivial term arises already at order
$\mathcal{O}\bigl{(}\Xi^{-2}\bigr{)}$ in the second part of Eq. (20) that
imposes the condition of fixed displacement, stemming from the corrections at
order $\mathcal{O}\bigl{(}\ell_{0}^{-2}\bigr{)}$ in Eqs. (16). Hence, already
at this order that does not even resolve the effect of non-zero differential
tension $\delta$, the governing equation differs from Euler’s Elastica
equation.. We therefore truncate the expansion at fourth order to obtain a
description of the epithelium in the spirit of a Landau theory. Not only do
nonlinear elastic terms arise at this order of truncation, but nonlocal terms
appear, too: the theory is _not_ elastic Libai and Simmonds (1998), since the
energy depends not only on strain (i.e. curvature), but also on its (spatial)
derivatives, introducing a nonlocal dependence on strain.

To obtain the governing equation, we vary the truncated expansion (20) with
respect to $\psi$, noting that $\Lambda$ is a constant since the trapezoidal
cells are required to match up exactly 444In Ref. Krajnc and Ziherl (2015),
the energy was expressed in terms of $\psi$ and $\ell$, which were then both
varied, while they are in fact related by $\ell\sec{\phi(\psi)}=\text{const.}$
for matching cells. If the lateral sides are not required to match up exactly,
corrections to the final term in Eq. (2) would have to be introduced, however,
to ensure a consistent description of the adhesion between neighboring cells.
Notwithstanding this, minimising Eq. (6) with respect to $\lambda$ leads to
$e(\phi)=2\ell_{0}\sqrt{\sec{\phi}+\delta\tan{\phi}}$, with
$e^{\prime\prime}(0)=\bigl{(}2-\delta^{2}\bigr{)}\ell_{0}/4$. For this reason,
if the lateral sides of the cells need not match up, the flat epithelium is
unstable to small perturbations if $\smash{\delta>\sqrt{2}}$, as obtained in
Ref. Krajnc and Ziherl (2015). After a considerable amount of algebra, we find

$\displaystyle\ddddot{\psi}=6\Xi^{2}\ddot{\psi}-\dfrac{3\Delta\Xi\dot{\psi}\ddot{\psi}}{1+\mu\cos{\psi}}+\dfrac{15+27\mu\cos{\psi}}{4(1+\mu\cos{\psi})}\dot{\psi}^{2}\ddot{\psi}+\dfrac{\mu\Xi^{4}\sin{\psi}}{1+\mu\cos{\psi}}\left[24-\dfrac{3}{\Xi^{2}}\dot{\psi}^{2}-\dfrac{1}{\Xi^{4}}\left(\dfrac{23}{16}\dot{\psi}^{4}-\dfrac{3}{2}\ddot{\psi}^{2}-2\dot{\psi}\dddot{\psi}\right)\right],$
(21)

wherein $\Delta=\delta\Lambda^{2}$, subject to the boundary conditions

$\displaystyle\psi(0)=\psi(1)=0,$
$\displaystyle\ddot{\psi}(0)=\ddot{\psi}(1)=0,$ (22a) and the integral
condition
$\displaystyle\int_{0}^{1}{\cos{\psi}\left(1+\dfrac{\dot{\psi}^{2}}{8\Xi^{2}}+\dfrac{41\dot{\psi}^{4}}{1920\Xi^{4}}+\dfrac{\ddot{\psi}^{2}}{40\Xi^{4}}+\dfrac{\dot{\psi}\dddot{\psi}}{240\Xi^{4}}\right)\mathrm{d}\sigma}=\Lambda(1-D).$
(22b)

The last condition imposes the fixed end-to-end shortening of the epithelium.
These equations have a trivial solution $\psi=0$, $\Lambda=(1-D)^{-1}$,
corresponding to the compressed but unbuckled state of the epithelium.

We note that, although Eq. (21) only depends on $\delta$ and $\Lambda$ through
their agglomerate $\Delta$, a separate dependence on $\Lambda$ arises in
condition (22b). Minimising the energy of buckled solutions of Eqs. (21,22)
with respect to $\Lambda$ finally determines $\Lambda$.

##  III Buckling Analysis

In this section, we analyze the buckling behavior of the epithelium, first
determining the threshold for buckling analytically and then discussing the
post-buckling behavior using a weakly nonlinear analysis of the governing
equations.

The buckling analysis naturally divides into two parts: we first seek buckled
configurations of small amplitude for each value of $\Lambda$, and then
minimize the energy of these configurations with respect to $\Lambda$.

###  III.1 Solution of the buckling problem

The form of the trivial solution and of condition (22b) suggest that the
appropriate small parameter for the first part of the analysis is

$\displaystyle\varepsilon^{2}=1-\Lambda(1-D).$ (23)

We therefore expand

$\displaystyle\psi(\sigma)=\varepsilon\Bigl{(}\psi_{0}(\sigma)+\varepsilon\psi_{1}(\sigma)+\varepsilon^{2}\psi_{2}(\sigma)+\mathcal{O}\bigl{(}\varepsilon^{3}\bigr{)}\Bigr{)},$
(24)

and write
$\mu=\mu_{0}+\varepsilon\mu_{1}+\varepsilon^{2}\mu_{2}+\mathcal{O}\bigl{(}\varepsilon^{3}\bigr{)}$.
It is important to note that, while the governing equations derived above are
only valid in the limit $\Xi\gg 1$, this parameter is not an asymptotic
parameter for the buckling analysis. It will however be useful to introduce

$\displaystyle\xi=\dfrac{\uppi}{\Xi}\ll 1.$ (25)

Next, we solve Eq. (21), subject to the boundary and integral conditions (22),
order-by-order.

####  III.1.1 Solution at order $\mathcal{O}(\varepsilon)$

At lowest order, the problem becomes

$\displaystyle\ddddot{\psi}_{0}-6\Xi^{2}\ddot{\psi}_{0}-\dfrac{24\Xi^{4}\mu_{0}}{1+\mu_{0}}\psi_{0}=0,$
(26)

subject to $\psi_{0}(0)=\psi_{0}(1)=0$,
$\ddot{\psi}_{0}(0)=\ddot{\psi}_{0}(1)=0$. The lowest eigenvalue of this
problem is

$\displaystyle\mu_{0}=\dfrac{z}{1-z},\quad\text{where
}z=\dfrac{\xi^{2}}{4}\left(1+\dfrac{\xi^{2}}{6}\right),$ (27)

the corresponding solution for $\psi_{0}$ being

$\displaystyle\psi_{0}(\sigma)=\Psi_{0}\sin{\uppi\sigma}.$ (28)

####  III.1.2 Solution at order $\mathcal{O}\bigl{(}\varepsilon^{2}\bigr{)}$

At next order, upon substituting for $\mu_{0}$,

$\displaystyle\dddot{\psi}_{1}-6\Xi^{2}\ddot{\psi}_{1}-24z\Xi^{4}\psi_{1}$
$\displaystyle\hskip
28.45274pt=24\mu_{1}\Xi^{4}(1-z)^{2}\psi_{0}-3\Delta\Xi(1-z)\dot{\psi}_{0}\ddot{\psi}_{0}.$
(29)

subject to $\psi_{1}(0)=\psi_{1}(1)=0$,
$\ddot{\psi}_{1}(0)=\ddot{\psi}_{1}(1)=0$. These conditions imply that
$\mu_{1}=0$, which is the usual result for the supercritical pitchfork
bifurcation expected for this buckling problem. Thence

$\displaystyle\psi_{1}(\sigma)=\Psi_{1}\sin{\uppi\sigma}+\Delta\Psi_{0}^{2}\left(\dfrac{\xi}{12}-\dfrac{13\xi^{3}}{144}\right)\sin{2\uppi\sigma}+\mathcal{O}\bigl{(}\xi^{5}\bigr{)}.$
(30)

####  III.1.3 Solution at order $\mathcal{O}\bigl{(}\varepsilon^{3}\bigr{)}$

Finally, upon substituting for $\mu_{0}$ and $\mu_{1}=0$,

$\displaystyle\ddddot{\psi}_{2}-6\Xi^{2}\ddot{\psi}_{2}-24z\Xi^{4}\psi_{2}$
$\displaystyle=24\Xi^{4}\mu_{2}(1-z)^{2}\psi_{0}-3\Delta\Xi(1-z)\Bigl{(}\dot{\psi}_{0}\ddot{\psi}_{1}+\ddot{\psi}_{0}\dot{\psi}_{1}\Bigr{)}+\left(\tfrac{15}{4}+3z\right)\dot{\psi}_{0}^{2}\ddot{\psi}_{0}$
$\displaystyle\hskip
113.81102pt-4z(1-3z)\Xi^{4}\psi_{0}^{3}-z\,\psi_{0}\left(3\Xi^{2}\dot{\psi}_{0}^{2}-\tfrac{3}{2}\ddot{\psi}_{0}^{2}-2\dot{\psi}_{0}\dddot{\psi}_{0}\right),$
(31)

subject to $\psi_{2}(0)=\psi_{2}(1)=0$,
$\ddot{\psi}_{2}(0)=\ddot{\psi}_{2}(1)=0$. After a considerable amount of
algebra, we obtain

$\displaystyle\mu_{2}=\Psi_{0}^{2}\left[\dfrac{\xi^{2}}{32}+\left(\dfrac{17}{384}-\dfrac{\Delta^{2}}{96}\right)\xi^{4}\right]+\mathcal{O}\bigl{(}\xi^{5}\bigr{)}.$
(32)

and thence

$\displaystyle\psi_{2}(\sigma)=\Psi_{2}\sin{\uppi\sigma}+\Psi_{0}\Psi_{1}\Delta\left(\dfrac{\xi}{6}-\dfrac{13\xi^{3}}{72}\right)\sin{2\uppi\sigma}+\Psi_{0}^{3}\left(\dfrac{1}{192}+\dfrac{4\Delta^{2}-9}{256}\xi^{2}-\dfrac{36\Delta^{2}-37}{768}\xi^{4}\right)\sin{3\uppi\sigma}+\mathcal{O}\bigl{(}\xi^{5}\bigr{)}.$
(33)

####  III.1.4 Calculation of the amplitudes

The constants $\Psi_{0},\Psi_{1},\Psi_{2}$ left undetermined by the above
calculation are obtained by expanding both sides of the integral condition
(22b). Solving order-by-order, we obtain

$\displaystyle\Psi_{0}$ $\displaystyle=\dfrac{2}{\sqrt{1-z}},$
$\displaystyle\Psi_{1}=0,$ (34)

where we have chosen $\Psi_{0}>0$ without loss of generality. The result
$\Psi_{1}=0$ is to be expected for a supercritical bifurcation; in fact, in a
standard elastic buckling problem, one would have $\psi_{1}\equiv 0$; here, a
non-zero $\psi_{1}$ is required because of the symmetry breaking resulting
from the term proportional to $\dot{\psi}^{3}$ in the Lagrangian (20).
Further, we obtain

$\displaystyle\Psi_{2}=\dfrac{1}{4}+\left(\dfrac{5}{32}-\dfrac{\Delta^{2}}{36}\right)\xi^{2}+\left(\dfrac{149}{7680}+\dfrac{37\Delta^{2}}{864}\right)\xi^{4}+\mathcal{O}\bigl{(}\xi^{5}\bigr{)}.$
(35)

###  III.2 Minimization of the Energy

In the second part of the buckling analysis, we determine the buckling
threshold, and then analyze the post-buckling behavior. Substituting for
$\Lambda$ using Eq. (23) in the energy term in Eq. (20) and expanding, the
energy of the buckled configuration is

$\displaystyle\left(\dfrac{1}{1-D}+1-D\right)+\left(\dfrac{1-D}{1-z}-\dfrac{1}{1-D}\right)\varepsilon^{2}+\mathcal{O}\bigl{(}\varepsilon^{3}\bigr{)},$
(36)

wherein the first bracketed term is the energy of the trivial solution
$\psi=0$, $\Lambda=(1-D)^{-1}$. Accordingly, buckled configurations become
energetically favorable if

$\displaystyle\dfrac{1-D}{1-z}-\dfrac{1}{1-D}<0\quad\Longleftrightarrow\quad
D>D_{\ast}\equiv 1-\sqrt{1-z}.$ (37)

In particular, the buckling threshold is independent of the differential
tension $\delta$.

We are left to determine the value of $\Lambda$ that minimizes the energy of
the buckled configuration. This is equivalent with relating $\varepsilon$ to
the excess compression $d=D-D_{\ast}>0$. We therefore write
$\varepsilon=\varepsilon_{0}d^{1/2}+\mathcal{O}(d)$, and obtain an expansion
of the energy in $d\ll 1$,

$\displaystyle\dfrac{\mathcal{E}}{\ell_{0}^{2}}=\left(C\varepsilon_{0}^{4}-\dfrac{2\varepsilon_{0}^{2}}{1-z}\right)d^{2}+\mathcal{O}\bigl{(}d^{5/2}\bigr{)},$
(38)

where

$\displaystyle
C=1+\dfrac{3\xi^{2}}{16}+\left(\dfrac{11}{64}-\dfrac{\delta^{2}}{48}\right)\xi^{4}+\mathcal{O}\bigl{(}\xi^{5}\bigr{)}.$
(39)

The energy is thus minimized when
$\varepsilon_{0}=\bigl{(}C(1-z)\bigr{)}^{-1/2}$. Substituting this result into
the expression for $\mu_{2}$, we finally obtain

$\displaystyle\mu\sim\mu_{0}(\xi)+\left(\dfrac{\xi^{2}}{8}+\dfrac{83-16\delta^{2}}{384}\xi^{4}+\mathcal{O}\bigl{(}\xi^{5}\bigr{)}\right)(D-D_{\ast}).$
(40)

This is the main result of our asymptotic analysis of the buckling: the force
required to compress the epithelium decreases with increasing differential
tension $\delta>0$.

In general, the buckled configuration features both energetically favorable
regions ($\dot{\psi}>0$ if $\delta>0$) and unfavorable regions
($\dot{\psi}<0$). Still, this result shows that buckling overall is
facilitated if $\delta>0$ compared to the $\delta=0$ case.

Figure 3: Numerical buckling results. (a) Plot of $\Lambda$ against relative
compression $D$. Above a critical compression $D_{\triangle}$ (buckled shape
for $D=D_{\triangle}$ shown), solution shapes at the energy minimum begin to
self-intersect (dotted line for $D>D_{\triangle}$). Inset: zoomed plot of
$\Lambda$ (filled marks) against $D$ close to the buckling threshold
$D_{\ast}$. (b) Plot of compressive force $\mu$ against relative compression
$D$, showing numerical results (solid line and dotted line for
$D>D_{\triangle}$) in agreement with asymptotic results (dashed line). Inset:
zoomed plot of $\mu$ against $D$ close to $D_{\ast}$. Parameter values for
numerical calculations: $\Xi=20$, $\delta=1$, $\ell_{0}=\sqrt{10}$. (c)
Critical compression $D_{\triangle}$ against $\smash{r=\Xi/\ell_{0}^{2}}$, for
different values of $\delta$, at fixed $\ell_{0}=\sqrt{10}$, and approximation
(44) thereof. Insets show buckled shapes at $D=D_{\triangle}$ and $\delta=1$,
for different values of $\Xi$.

##  IV Post-buckling behavior

While asymptotic analysis can describe the deformations of the epithelium just
beyond the buckling threshold, larger compressions must be studied
numerically. We solve the governing equation (21), complemented by the
boundary and integral conditions (22), numerically using the boundary-value-
problem solvers bvp4c, bvp5c of Matlab (The MathWorks, Inc.) and the
continuation package auto Doedel  _et al._ (2012).

###  IV.1 Transition to Constricted Cells

For the numerical solution, we fix $\Xi$ and $\delta$, and obtain solutions
for different values of $\Lambda$. By interpolation, we determine the value of
$\Lambda$ that minimizes the energy (Fig. 3a). Thence we obtain the
corresponding value of the compressive force $\mu$ (Fig. 3b), in agreement
with the asymptotic results of the previous section. This also validates our
numerical implementation of the system.

There is, however, one extra constraint that has not been incorporated into
the continuum equations: the constraint, related to the transition to
constricted cells that we have briefly discussed before, that the lateral
sides of the trapezoidal cells cannot self-intersect. At the level of the
continuum description, this constraint translates to the condition that the
apical and basal surfaces of the epithelium cannot self-intersect. The apical
and basal surfaces of the cell sheet are described by the curves

$\displaystyle x_{\pm}$ $\displaystyle=x\mp\ell\sin{\psi},$ $\displaystyle
y_{\pm}$ $\displaystyle=y\pm\ell\cos{\psi},$ (41)

where $\ell=\ell_{0}\Lambda\cos{\phi}$ is the local thickness of the cell
sheet. It is important to note that Eqs. (41) are exact equations since the
Kirchhoff ‘hypothesis’ of the analysis of slender elastic structures, the
asymptotic result Audoly and Pomeau (2010) that the normal to the undeformed
midline remains normal in the deformed configuration, is an exact result in
the discrete model that underlies our continuum theory. (For this same reason,
an analogous analysis for an elastic object beyond asymptotically small
deformations would, rather more intricately, require solving for the stretches
in each parallel to the midline.)

Numerically, we find that, as $D$ is increased at fixed $\Lambda$, the shapes
of minimal energy self-intersect above a critical compression $D_{\triangle}$.
The numerical solutions also reveal that self-intersections first arise at
$\sigma=0$, when $\dot{x}_{+}=0$ there. Expanding this condition using
$\psi(0)=0$ and $\dot{\phi}(0)=0$, of which the latter follows from Eq. (13)
by symmetry,

$\displaystyle\dfrac{f(0)}{\Lambda}-\Lambda\dfrac{\ell_{0}^{2}}{\Xi}\dot{\psi}(0)\cos{\phi(0)}=0,$
(42)

where $f$ is defined as in Eq. (16a). We note that, with this condition, an
explicit dependence on both $\ell_{0}$ and $\Xi$ has arisen for the first time
in our analysis.

Figure 4: Buckling for $D>D_{\triangle}$. Two scenarios are possible: (a) If
$\partial\dot{x}(0)/\partial\Lambda>0$ or $\Xi<\Xi_{\ast}(\delta)$, buckled
shapes with increased $\Lambda$ do not self-intersect. (b) If
$\partial\dot{x}(0)/\partial\Lambda<0$ or $\Xi>\Xi_{\ast}(\delta)$, buckled
shapes with decreased $\Lambda$ do not self-intersect. Numerical results: (c)
Plot of $\Lambda$ (for buckled shapes of minimal energy without self-
intersections) against relative compression $D$, for $\Xi_{1}<\Xi_{\ast}$
(thick lines) and $\Xi_{2}>\Xi_{\ast}$ (thin lines). For $\Xi<\Xi_{\ast}$,
$\Lambda\rightarrow(1-D)^{-1}$ for $D>D_{\triangle}$. Insets show buckled
shapes at $D=D_{\triangle}$. (d) Corresponding plots of scaled compressive
force $\mu/\xi^{2}$ against $D$. Dotted lines for $D>D_{\triangle}$ correspond
to self-intersecting shapes at the energy minimum. Parameter values for
numerical calculations: $\delta=1$, $\Xi_{1}=20$, $\Xi_{2}=34$,
$\ell_{0}=\sqrt{10}$.

#### Estimating the critical compression $D_{\triangle}$

The numerical data in Fig. 3a,b suggest that the asymptotic results of the
previous section can approximate the buckling behavior of the epithelium well
up to compressions as large as $D_{\triangle}$. We therefore use our
asymptotic results to estimate the critical compression $D_{\triangle}$. For
this purpose, we treat $r=\smash{\Xi/\ell_{0}^{2}}$ as an $\mathcal{O}(1)$
quantity. Then, using $\Lambda=(1-z)^{-1/2}+\mathcal{O}(d)$,

$\displaystyle\dfrac{f(0)}{\Lambda}-\Lambda\dfrac{\ell_{0}^{2}}{\Xi}\dot{\psi}(0)\cos{\phi(0)}=\sqrt{1-z}-\dfrac{\uppi\varepsilon_{0}\Psi_{0}}{r\sqrt{1-z}}d^{1/2}+\mathcal{O}(d),$
(43)

whence, to lowest order in $\xi$,

$\displaystyle D_{\triangle}\approx\dfrac{r^{2}}{4\uppi^{2}}.$ (44)

This approximation is not itself an asymptotic result, yet, for small enough
values of $r$, it compares well to numerical estimates of $D_{\triangle}$
obtained by a bisection search (Fig. 3c). The numerical results also show
that, at fixed $\Xi$, $D_{\triangle}$ decreases with increasing $\delta$, with
relative variations of about $10\%$ for the range of $\delta$ under
consideration.

We also find numerically that, at large values of $r$, steric interactions
between different parts of the cell sheet become important before $D$ reaches
$D_{\triangle}$. A detailed analysis of these interactions is beyond the scope
of this discussion.

We have tacitly assumed that, at fixed $\Xi$ and at fixed $D>D_{\ast}$, the
energy $\mathcal{E}$ has a single local minimum as a function of $\Lambda$.
This is indeed the case for small enough values of $\delta$, but fails at
compressions $D>D_{\triangle}$ as $\delta$ is increased, so this possibility
is not of direct relevance to the present discussion. Interestingly,
eigenmodes (buckled solutions with zero force) of the epithelium arise at
large $\delta$. These eigenmodes are not energy minimizers, but, for
completeness, we discuss these solutions in Appendix B.

###  IV.2 Buckled shapes for $\boldsymbol{D>D_{\triangle}}$

As $D$ is increased beyond $D_{\triangle}$, we might expect fans of
constricted cells to expand around the trough and (later) the crest of the
buckled shape, but deriving the equations describing these fans and solving
for these shapes is beyond the scope of the present discussion.

Here, we note simply that, for $D>D_{\triangle}$, buckled shapes without self-
intersections can be found. Two scenarios can be envisaged _a priori_ ,
depending on the sign of $\partial\dot{x}(0)/\partial\Lambda$ at the energy
minimum (Fig. 4a,b): if $\partial\dot{x}(0)/\partial\Lambda>0$, buckled
solutions without self-intersection arise as $\Lambda$ is increased; if
$\partial\dot{x}(0)/\partial\Lambda>0$, such shapes are found as $\Lambda$ is
decreased. Interestingly, both of these possibilities arise in the system:
there exists a critical value $\Xi=\Xi_{\ast}(\delta)$ at which
$\partial\dot{x}(0)/\partial\Lambda=0$. The first possibility occurs in the
case $\Xi<\Xi_{\ast}$, and the second one in the case $\Xi>\Xi_{\ast}$ (Fig.
4c).

If $\Xi<\Xi_{\ast}$, the buckling amplitude decreases for these solutions as
$D$ is increased beyond $D_{\triangle}$: $\Lambda$ tends to its value
$(1-D)^{-1}$ in the unbuckled, compressed configuration, while $\mu$ decreases
(Fig. 4c,d). By contrast, if $\Xi>\Xi_{\ast}$, the buckling amplitude is
increased: as $D$ increases and $\Lambda$ decreases, $\mu$ increases faster
than in the self-intersecting configurations at the energy minimum (Fig.
4c,d).

While these qualitative considerations cannot capture the exact mechanics of
the fans of constricted cells near the through and crest of the buckled shape,
we expect them to give a qualitative indication of the buckling behavior as
$D$ is increased just beyond $D_{\triangle}$. For larger values of $D$, there
are more intricate possibilities: there are in general two values of $\Lambda$
such that $\dot{x}(0)=0$, on either side of the energy minimum. One of these
solutions defines the branch shown in Fig. 4c, but the second solution may
become energetically favorable over the first one as $D$ is increased
sufficiently. This actually happens on the branch with
$\Xi=\Xi_{1}<\Xi_{\ast}$ in Fig. 4c at $D\approx 0.44$, but, for the second
solution, different parts of the cell sheet start to touch before this value
of $D$ is reached and hence we do not pursue this further here.

##  V Conclusion

In this paper, we have derived, by taking a rigorous asymptotic limit, the
continuum limit of a simple discrete differential-tension model of a two-
dimensional epithelium. If the expansion is carried to high enough order for
the differential tension between the apical and basal sides of the epithelium
to arise in the energy, nonelastic terms that are nonlocal in the strains
appear. This is the key lesson to be drawn from taking the continuum limit. We
have gone on to use this continuum model to study the buckling of the
epithelium under imposed confinement, showing how, post-buckling, the
compressive force is reduced with increasing differential tension. A second
buckling transition occurs when constricted cells start to form near the
troughs and crests of the buckled shape; we have discussed the behavior close
to this transition qualitatively. Taking the analysis of the buckling behavior
of epithelium in this continuum framework beyond the transition to constricted
cells is the key challenge for future work on this problem.

Possible extensions of the continuum framework include mimicking the setup of
studies of the discrete model Štorgel  _et al._ (2016); Brezavšček  _et al._
(2012); Rauzi  _et al._ (2015) by coupling the epithelium to an elastic
substrate or incorporating fixed-volume constraints for a closed one-
dimensional epithelium. The question of how to extend the continuum model to
describe a two-dimensional epithelium also remains open. In particular, how
are the deviations from elasticity affected by the increase of the
dimensionality of the system?

Cell sheet deformations during development commonly feature large _geometric_
deformations, but the _elastic_ deformations can remain small provided that
the deformed geometry remains close to the intrinsic geometry that is
generally different from the initial geometry because of cell shape changes,
cell intercalation, and related processes. For this reason, developmental
events as intricate as the inversion process of the green alga _Volvox_ can be
modelled quantitatively using a Hookean shell theory Höhn  _et al._ (2015);
Haas  _et al._ (2018). By contrast, large deformations of many a biological
material are not in general described well by neo-Hookean constitutive
equations, although other families of hyperelastic constitutive equations
predict behavior in quantitative agreement with experimental data for brain
and fat tissues Mihai  _et al._ (2015, 2017); Budday  _et al._ (2017).
However, and in spite of the ubiquity of these elastic models, in particular
in the modelling of the folds of the cerebrum Richman  _et al._ (1975);
Tallinen  _et al._ (2013); Manyuhina  _et al._ (2014); Budday  _et al._
(2015); Tallinen  _et al._ (2016), it was pointed out very recently that the
folding of the cerebellum is fundamentally inconsistent with the differential-
growth hypothesis Engstrom  _et al._ (2018): in the cerebellum, the
oscillations of the thicknesses of the core and the growing cortex are out-of-
phase, while elastic bilayer instabilities lead to in-phase oscillations
Engstrom  _et al._ (2018). All of this emphasizes the need for a deeper
understanding of how continuum models relate to properties of structures at
the cell level. By explicitly showing how both nonlinear and nonlocal elastic
terms arise in the continuum limit of a simple discrete model and impact on
its behavior, the present analysis has taken a first step in this direction.

###### Acknowledgements.

We thank Matej Krajnc for discussions of Ref. Krajnc and Ziherl (2015) at an
early stage of this work, and are grateful for support from the Engineering
and Physical Sciences Research Council (Established Career Fellowship
EP/M017982/1, REG; Doctoral Prize Fellowship, PAH), Wellcome Trust
Investigator Award 207510/Z/17/Z, the Schlumberger Chair Fund, and Magdalene
College, Cambridge (Nevile Research Fellowship, PAH).

##  Appendix A Derivation of Eq. (12)

In this appendix, we relate $\phi$ and $\psi$ by deriving Eq. (12) used in the
main text. At the same time, we verify that this expansion is indeed
consistent at all orders. Taylor expanding the right-hand side of Eq. (11),

$\displaystyle\psi(s+ks_{0})-\psi(s)$
$\displaystyle=\phi(s)+2\left(\sum_{j=1}^{k-1}{\phi(s+js_{0})}\right)+\phi(s+ks_{0})=2k\phi(s)+\sum_{n=1}^{\infty}{\dfrac{s_{0}^{n}}{n!}\phi^{(n)}(s)\left\\{k^{n}+2\sum_{j=1}^{k-1}{j^{n}}\right\\}}$
$\displaystyle\hskip
71.13188pt=2k\phi(s)+\sum_{n=1}^{\infty}{\dfrac{s_{0}^{n}}{n!}\phi^{(n)}(s)\left\\{k^{n}+\dfrac{2}{n+1}\sum_{j=0}^{n}{(-1)^{j}\binom{n+1}{j}\mathcal{B}_{j}(k-1)^{n+1-j}}\right\\}},$
(A1)

where we have used Faulhaber’s formula Conway and Guy (1996) to expand the sum
of powers of integers, and where
$\mathcal{B}_{0}=1,\mathcal{B}_{1}=-\frac{1}{2},\dots$ denote the Bernoulli
numbers (of the first kind) Spanier and Oldham (1987). Expanding
$\smash{(k-1)^{n+1-j}}$ using the binomial theorem and simplifying the
binomial coefficients,

$\displaystyle\psi(s+ks_{0})-\psi(s)=2k\phi(s)+\sum_{n=1}^{\infty}{\dfrac{\phi^{(n)}(s)}{n!\ell_{0}^{n}}A(k,n)},$
(A2)

where we have introduced

$\displaystyle A(k,n)$
$\displaystyle=k^{n}+2(-1)^{n+1}n!\sum_{j=0}^{n}{\sum_{i=0}^{n+1-j}{\dfrac{(-1)^{i}\mathcal{B}_{j}k^{i}}{i!j!(n+1-i-j)!}}}=\sum_{i=0}^{n+1}{a_{i}(n)k^{i}},$
(A3)

wherein the coefficients $a_{0},a_{1},\dots,a_{n+1}$ depend on $n$. We notice
in particular that

$\displaystyle a_{0}$ $\displaystyle=0,$ $\displaystyle a_{n}=0,$
$\displaystyle a_{n+1}=\dfrac{2}{n+1},$ (A4)

of which the last two are obtained by direct computation, and the first one
follows using an identity of Bernoulli numbers Spanier and Oldham (1987),

$\displaystyle\sum_{j=0}^{n}{\binom{n+1}{j}\mathcal{B}_{j}}=0\qquad\mbox{for
}n=1,2,\dots.$ (A5)

Moreover, for $i=1,2,\dots,n-1$,

$\displaystyle
a_{i}(n)=2(-1)^{n+1-i}\dfrac{n!}{i!}\sum_{j=0}^{n+1-i}{\dfrac{\mathcal{B}_{j}}{j!(n+1-i-j)!}}.$
(A6)

Accordingly,

$\displaystyle
A(k,n)=\dfrac{2k^{n+1}}{n+1}+2(-1)^{n+1}n!\sum_{i=1}^{n-1}{(-1)^{i}\dfrac{k^{i}}{i!}\left\\{\sum_{j=0}^{n+1-i}{\dfrac{\mathcal{B}_{j}}{j!(n+1-i-j)!}}\right\\}}.$
(A7)

Upon inverting the order of summation, Eq. (A2) becomes

$\displaystyle\psi(s+ks_{0})-\psi(s)=2k\phi(s)+\sum_{n=1}^{\infty}{\dfrac{\phi^{(n)}(s)}{n!\ell_{0}^{n}}\left(\dfrac{2k^{n+1}}{n+1}\right)}+\sum_{i=1}^{\infty}{\sum_{n=i+1}^{\infty}{2(-1)^{n+1-i}\dfrac{\phi^{(n)}(s)}{\ell_{0}^{n}}}\dfrac{k^{i}}{i!}\left\\{\sum_{j=0}^{n+1-i}{\dfrac{\mathcal{B}_{j}}{j!(n+1-i-j)!}}\right\\}}.$
(A8)

Finally, upon relabelling indices in the first summation and changing
variables $n\longmapsto m=n+1-i$ in the last summation,

$\displaystyle\psi(s+ks_{0})-\psi(s)=\sum_{i=1}^{\infty}{\dfrac{k^{i}}{i!\ell_{0}^{i}}\left\\{2\ell_{0}\phi^{(i-1)}(s)+\sum_{m=2}^{\infty}{\dfrac{2(-1)^{m}}{\ell_{0}^{m-1}}}\phi^{(i-1+m)}(s)\left(\sum_{j=0}^{m}{\dfrac{\mathcal{B}_{j}}{j!(m-j)!}}\right)\right\\}}.$
(A9)

But, rearranging Eq. (A5),

$\displaystyle\sum_{j=0}^{m}{\dfrac{\mathcal{B}_{j}}{j!(m-j)!}}=\dfrac{\mathcal{B}_{m}}{m!}\qquad\mbox{for
}m=2,3,\dots.$ (A10)

Since $\mathcal{B}_{n}=0$ for odd $n>1$, and using $\mathcal{B}_{0}=1$, we
finally obtain

$\displaystyle\psi(s+ks_{0})-\psi(s)=\sum_{i=1}^{\infty}{\dfrac{k^{i}}{i!\ell_{0}^{i}}\left\\{\sum_{m=0}^{\infty}{\dfrac{2\mathcal{B}_{2m}}{(2m)!}\dfrac{\phi^{(i-1+2m)}(s)}{\ell_{0}^{2m-1}}}\right\\}}.$
(A11)

Comparing this to the Taylor expansion of the left-hand side,

$\displaystyle\psi(s+ks_{0})-\psi(s)=\sum_{i=1}^{\infty}{\dfrac{k^{i}}{i!\ell_{0}^{i}}\psi^{(i)}(s)},$
(A12)

we deduce that the expansion is consistent at all orders, with, in particular,

$\displaystyle\psi^{\prime}(s)=\sum_{m=0}^{\infty}{\psi_{m}\dfrac{\phi^{(2m)}(s)}{\ell_{0}^{2m-1}}}\qquad\text{where
}\psi_{m}=\dfrac{2\mathcal{B}_{2m}}{(2m)!},$ (A13)

which is Eq. (12). As noted in the main text, we are not aware of a closed
form for the inverted series that expresses $\phi$ as a function of the
derivatives of $\psi$. Formally, inverting Eq. (12) gives

$\displaystyle\phi(s)=\sum_{m=0}^{\infty}{\phi_{m}\dfrac{\psi^{(2m+1)}(s)}{\ell_{0}^{2m+1}}},$
(A14)

where the coefficients $\phi_{0},\phi_{1},\dots$ are determined recursively by
$\phi_{0}\psi_{0}=1$ and

$\displaystyle\sum_{j=0}^{m}{\phi_{j}\psi_{m-j}}=0\qquad\text{for
}m=1,2,\dots.$ (A15)

In agreement with Eq. (13), we find

$\displaystyle\phi_{0}=\tfrac{1}{2},$ $\displaystyle\phi_{1}=-\tfrac{1}{24},$
$\displaystyle\phi_{2}=\tfrac{1}{240},$ $\displaystyle\dots.$ (A16)

Figure 5: Eigenmodes of a buckled epithelium: plot of relative end-to-end
shortening $D$ against $\Delta$. Parameter value: $\Xi=20$. No eigenmodes were
found for $\Delta<\Delta_{\ast}(\Xi)$. On the dashed part of the branch,
$\mathcal{E}<2$. Continuation failed at the point marked $\times$. Inset: plot
of $\Delta_{\ast}$ against $\Xi$ (solid line) and power-law fit (dashed line).

##  Appendix B Eigenmodes of the buckled epithelium

Eigenmodes of the epithelium are non-zero solutions of the governing equation
(21) with $\mu=0$. They thus obey

$\displaystyle\ddddot{\psi}=6\Xi^{2}\ddot{\psi}-3\Delta\Xi\dot{\psi}\ddot{\psi}+\dfrac{15}{4}\dot{\psi}^{2}\ddot{\psi},$
(B1)

subject to

$\displaystyle\psi(0)=\psi(1)=0,$
$\displaystyle\ddot{\psi}(0)=\ddot{\psi}(1)=0.$ (B2)

To find eigenmodes numerically, we remove the trivial, zero solution by
imposing a non-zero compression $D$ and varying this compression until a
solution with $\mu=0$ is found.

Numerically, we obtain eigenmodes if $\Delta\geqslant\Delta_{\ast}$, but find
no solutions if $\Delta<\Delta_{\ast}$, for some value $\Delta_{\ast}$
depending on $\Xi$ (Fig. 5). Some of these solutions have energy
$\mathcal{E}<2$ (Fig. 5), lower than the energy of the uncompressed, flat
solution; these are spontaneous buckled modes that arise in the absence of
external forces, but, as is apparent from the corresponding values $D>1$ (Fig.
5), these solutions are unphysical. Plotting $\Delta_{\ast}$ against $\Xi$
(Fig. 5, inset) suggests that $\Delta_{\ast}$ approaches a constant value as
$\Xi$ grows large. We observe that the numerical data are well approximated by
a power-law $\Delta_{\ast}=c_{1}+c_{2}\Xi^{-5/4}$, where $c_{1}\approx 3.96$,
$c_{2}\approx 19.5$ (Fig. 5, inset).

The ‘large’ values of $\Delta$ and hence $\delta$ for these eigenmodes beckon
a comment on the formal range of validity of the continuum model: stability of
the underlying discrete model requires $\alpha,\beta\geqslant 0$ Krajnc  _et
al._ (2013), and hence $\delta\leqslant\ell_{0}^{2}$. While the asymptotic
expansion leading to the geometric relation (13) was an expansion in the large
parameter $\ell_{0}$, it did not involve $\delta$. By contrast, the expansion
of the Lagrangian (20), which did involve $\delta$, was an expansion in a
different large parameter, $\Xi$. Hence large values of
$\delta\lesssim\ell_{0}^{2}$ are indeed in the formal range of validity of the
continuum model provided that $\Xi$ is large enough.

## References

  * Lecuit and Lenne (2007) T. Lecuit and P.-F. Lenne, “Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis,” Nat. Rev. Mol. Cell Biol. 8, 633–644 (2007).
  * Odell  _et al._ (1981) G.M. Odell, G. Oster, P. Alberch, and B. Burnside, “The mechanical basis of morphogenesis,” Dev. Biol. 85, 446–462 (1981).
  * Rauzi  _et al._ (2013) M. Rauzi, A. Hočevar Brezavšček, P. Ziherl, and M. Leptin, “Physical models of mesoderm invagination in _Drosophila_ embryo,” Biophys. J. 105, 3–10 (2013).
  * Fletcher  _et al._ (2016) A. G. Fletcher, F. Cooper, and R. E. Baker, “Mechanocellular models of epithelial morphogenesis,” Phil. Trans. R. Soc. B 372, 20150519 (2016).
  * Derganc  _et al._ (2009) J. Derganc, S. Svetina, and B. Žekš, “Equilibrium mechanics of monolayered epithelium,” J. Theor. Biol. 260, 333–339 (2009).
  * Brezavšček  _et al._ (2012) A. Hočevar Brezavšček, M. Rauzi, M. Leptin, and P. Ziherl, “A model of epithelial invagination driven by collective mechanics of identical cells,” Biophys. J. 103, 1069–1077 (2012).
  * Krajnc  _et al._ (2013) M. Krajnc, N. Štorgel, A. Hočevar Brezavšček, and P. Ziherl, “A tension-based model of flat and corrugated simple epithelia,” Soft Matter 9, 8368–8377 (2013).
  * Hannezo  _et al._ (2014) E. Hannezo, J. Prost, and J.-F. Joanny, “Theory of epithelial sheet morphology in three dimensions,” Proc. Nat. Acad. Sci. USA 111, 27–32 (2014).
  * Rauzi  _et al._ (2015) M. Rauzi, U. Krzic, T. E. Saunders, M. Krajnc, P. Ziherl, L. Hufnagel, and M. Leptin, “Embryo-scale tissue mechanics during _Drosophila_ gastrulation movements,” Nat. Commun. 6, 8677 (2015).
  * Krajnc and Ziherl (2015) M. Krajnc and P. Ziherl, “Theory of epithelial elasticity,” Phys. Rev. E 92, 052713 (2015).
  * Štorgel  _et al._ (2016) N. Štorgel, M. Krajnc, P. Mrak, J. Štrus, and P. Ziherl, “Quantitative morphology of epithelial folds,” Biophys. J. 110, 269–277 (2016).
  * Krajnc  _et al._ (2018) M. Krajnc, S. Dasgupta, P. Ziherl, and J. Prost, “Fluidization of epithelial sheets by active cell rearrangements,” Phys. Rev. E 98, 022409 (2018).
  * Howard  _et al._ (2011) J. Howard, S. W. Grill, and J. S. Blois, “Turing’s next steps: the mechanochemical basis of morphogenesis,” Nat. Rev. Mol. Cell Biol. 12, 392–398 (2011).
  * Dasbiswas  _et al._ (2018) K. Dasbiswas, E. Hannezo, and N. S. Gov, “Theory of epithelial cell shape transitions induced by mechanoactive chemical gradients,” Biophys. J. 114, 968–977 (2018).
  * Richman  _et al._ (1975) D. P. Richman, R. M. Stewart, J. W. Hutchinson, and V. S. Caviness, Jr., “Mechanical model of brain convolutional development,” Science 189, 18–21 (1975).
  * Tallinen  _et al._ (2013) T. Tallinen, J. S. Biggins, and L. Mahadevan, “Surface sulci in squeezed soft solids,” Phys. Rev. Lett. 110, 024302 (2013).
  * Manyuhina  _et al._ (2014) O. V. Manyuhina, D. Mayett, and J. M. Schwarz, “Elastic instabilities in a layered cerebral cortex: a revised axonal tension model for cortex folding,” New J. Phys. 16, 123058 (2014).
  * Budday  _et al._ (2015) S. Budday, P. Steinmann, and E. Kuhl, “Secondary instabilities modulate cortical complexity in the mammalian brain,” Philos. Mag. 95, 3244–3256 (2015).
  * Tallinen  _et al._ (2016) T. Tallinen, J. Y. Chung, F. Rousseau, N. Girard, J. Lefève, and L. Mahadevan, “On the growth and form of cortical convolutions,” Nat. Phys. 12, 588–593 (2016).
  * Lejeune  _et al._ (2016) E. Lejeune, A. Javili, J. Weickenmeier, E. Kuhl, and C. Linder, “Tri-layer wrinkling as a mechanism for anchoring center initiation in the developing cerebellum,” Soft Matter 12, 5613–5620 (2016).
  * Hannezo  _et al._ (2011) E. Hannezo, J. Prost, and J.-F. Joanny, “Instabilities of monolayered epithelia: Shape and structure of villi and crypts,” Phys. Rev. Lett. 107, 078104 (2011).
  * Savin  _et al._ (2011) T. Savin, N. A. Kurpios, A. E. Shyer, P. Florescu, H. Liang, L. Mahadevan, and C. J. Tabin, “On the growth and form of the gut,” Nature 476, 57–62 (2011).
  * Shyer  _et al._ (2013) A. E. Shyer, T. Tallinen, N. L. Nerurkar, Z. Wei, E. S. Gil, D. L. Kaplan, C. J. Tabin, and L. Mahadevan, “Villification: how the gut gets its villi,” Science 342, 212–218 (2013).
  * Hufnagel  _et al._ (2007) L. Hufnagel, A. A. Teleman, H. Rouault, S. M. Cohen, and B. I. Shraiman, “On the mechanism of wing size determination in fly development,” Proc. Nat. Acad. Sci. USA 104, 3835–3840 (2007).
  * Farhadifar  _et al._ (2007) R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, “The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing,” Curr. Biol. 17, 2095–2104 (2007).
  * Hilgenfeldt  _et al._ (2008) S. Hilgenfeldt, S. Erisken, and R. W. Carthew, “Physical modeling of cell geometric order in an epithelial tissue,” Proc. Nat. Acad. Sci. USA 105, 907–911 (2008).
  * Bi  _et al._ (2015) D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning, “A density-independent rigidity transition in biological tissues,” Nat. Phys. 11, 1074–1079 (2015).
  * Bi  _et al._ (2016) D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, “Motility-driven glass and jamming transitions in biological tissues,” Phys. Rev. X 6, 021011 (2016).
  * Landau and Lifshitz (1970) L. D. Landau and E. M. Lifshitz,  _Theory of Elasticity_ , 2nd ed., Course of Theoretical Physics, Vol. 7 (Pergamon, Oxford, England, 1970) Chap. 21, pp. 97–100.
  * Audoly and Pomeau (2010) B. Audoly and Y. Pomeau,  _Elasticity and Geometry_ (Oxford University Press, Oxford, England, 2010).
  * Note (1) In this respect, the present analysis differs from that of Ref. Krajnc and Ziherl (2015), where the energy was expressed as an integral with respect to the arclength $S$ in the deformed configuration, related in fact to $s$ by $\ell\tmspace+{.1667em}\mathrm{d}S=\ell_{0}\tmspace+{.1667em}\mathrm{d}s$. In Ref. Krajnc and Ziherl (2015), a weaker global area conservation constraint for the cell sheet had therefore to be imposed separately.
  * Spanier and Oldham (1987) J. Spanier and K. B. Oldham,  _An Atlas of Functions_ (Hemisphere, Washington, DC, USA, 1987) Chap. 4, pp. 35–38.
  * Note (2) These and other expansions in this paper were carried out with the help of Mathematica (Wolfram, Inc.) to facilitate the manipulation of complicated algebraic expressions.
  * Note (3) A non-trivial term arises already at order $\mathcal{O}\mathopen{\hbox{\frozen@everymath$\left(\vbox to1.0pt{}\right.$}}\Xi^{-2}\mathclose{\hbox{\frozen@everymath$\left)\vbox to1.0pt{}\right.$}}$ in the second part of Eq. (20) that imposes the condition of fixed displacement, stemming from the corrections at order $\mathcal{O}\mathopen{\hbox{\frozen@everymath$\left(\vbox to1.0pt{}\right.$}}\ell_{0}^{-2}\mathclose{\hbox{\frozen@everymath$\left)\vbox to1.0pt{}\right.$}}$ in Eqs. (16). Hence, already at this order that does not even resolve the effect of non-zero differential tension $\delta$, the governing equation differs from Euler’s Elastica equation.
  * Libai and Simmonds (1998) A. Libai and J. G. Simmonds,  _The Nonlinear Theory of Elastic Shells_, 2nd ed. (Cambridge University Press, Cambridge, England, 1998) Chap. III.J, pp. 36–38.
  * Note (4) In Ref. Krajnc and Ziherl (2015), the energy was expressed in terms of $\psi$ and $\ell$, which were then both varied, while they are in fact related by $\ell\mathop{sec}\nolimits{\phi(\psi)}=\text{const.}$ for matching cells. If the lateral sides are not required to match up exactly, corrections to the final term in Eq. (2) would have to be introduced, however, to ensure a consistent description of the adhesion between neighboring cells. Notwithstanding this, minimising Eq. (6) with respect to $\lambda$ leads to $e(\phi)=2\ell_{0}\sqrt{\mathop{sec}\nolimits{\phi}+\delta\mathop{tan}\nolimits{\phi}}$, with $e^{\prime\prime}(0)=\mathopen{\hbox{\frozen@everymath$\left(\vbox to1.0pt{}\right.$}}2-\delta^{2}\mathclose{\hbox{\frozen@everymath$\left)\vbox to1.0pt{}\right.$}}\ell_{0}/4$. For this reason, if the lateral sides of the cells need not match up, the flat epithelium is unstable to small perturbations if $\smash{\delta>\sqrt{2}}$, as obtained in Ref. Krajnc and Ziherl (2015).
  * Doedel  _et al._ (2012) E. J. Doedel, B. E. Oldman, A. R. Champneys, F. Dercole, T. Fairgrieve, Y. Kuznetsov, R. Paffenroth, B. Sandstede, X. Wang, and C. Zhang,  _Auto-07p : Continuation and Bifurcation Software for Ordinary Differential Equations_, Tech. Rep. (Concordia University, Montreal, Canada, 2012).
  * Höhn  _et al._ (2015) S. Höhn, A. R. Honerkamp-Smith, P. A. Haas, P. Khuc Trong, and R. E. Goldstein, “Dynamics of a _Volvox_ embryo turning itself inside out,” Phys. Rev. Lett. 114, 178101 (2015).
  * Haas  _et al._ (2018) P. A. Haas, S. S. M. H. Höhn, A. R. Honerkamp-Smith, J. B. Kirkegaard, and R. E. Goldstein, “The noisy basis of morphogenesis: mechanisms and mechanics of cell sheet folding inferred from developmental variability,” PLOS Biol. 16, e2005536 (2018).
  * Mihai  _et al._ (2015) L. A. Mihai, L. Chin, P. A. Janmey, and A. Goriely, “A comparison of hyperelastic constitutive models applicable to brain and fat tissues,” J. R. Soc. Interface 12, 20150486 (2015).
  * Mihai  _et al._ (2017) L. A. Mihai, S. Budday, G. A. Holzapfel, E. Kuhl, and A. Goriely, “A family of hyperelastic models for human brain tissue,” J. Mech. Phys. Solids 106, 60–79 (2017).
  * Budday  _et al._ (2017) S. Budday, G. Sommer, C. Birkl, C. Langkammer, J. Haybaeck, J. Kohnert, M. Bauer, F. Paulsen, P. Steinmann, E. Kuhl, and G. A. Holzapfel, “Mechanical characterization of human brain tissue,” Acta Biomater. 48, 319–340 (2017).
  * Engstrom  _et al._ (2018) T. A. Engstrom, T. Zhang, A. K. Lawton, A. L. Joyner, and J. M. Schwarz, “Buckling without bending: a new paradigm in morphogenesis,” (2018), arXiv:1806.06961 .
  * Conway and Guy (1996) J. H. Conway and R. K. Guy,  _The Book of Numbers_ (Springer, New York, NY, USA, 1996) p. 107.


# MLRM: A Multiple Linear Regression based Model for Average Temperature
Prediction of A Day

Ishu Gupta*
Cloud Computing Research Center
Department of Computer Science and Engineering
National Sun Yat-sen University
Kaohsiung, Taiwan
<EMAIL_ADDRESS>
Harsh Mittal
Department of Computer Applications
National Institute of Technology
Kurukshetra, India
136119
<EMAIL_ADDRESS>Deepak Rikhari
Department of Computer Applications
National Institute of Technology
Kurukshetra, India
136119
<EMAIL_ADDRESS>Ashutosh Kumar Singh
Department of Computer Applications
National Institute of Technology
Kurukshetra, India
136119
<EMAIL_ADDRESS>

###### Abstract

Weather is a phenomenon that affects everything and everyone around us on a
daily basis. Weather prediction has been an important point of study for
decades as researchers have tried to predict the weather and climatic changes
using traditional meteorological techniques. With the advent of modern
technologies and computing power, we can do so with the help of machine
learning techniques. We aim to predict the weather of an area using past
meteorological data and features using the Multiple Linear Regression Model.
The performance of the model is evaluated and a conclusion is drawn. The model
is successfully able to predict the average temperature of a day with an error
of 2.8 degrees Celsius.

_Keywords_ Weather forecasting $\cdot$ Machine learning $\cdot$ Multiple
linear regression $\cdot$ Prediction $\cdot$ Artificial neural network

##  1 Introduction

Weather prediction is the science of predicting the conditions of the
atmosphere at a particular time and place [1, 2, 3]. This prediction has been
attempted formally and informally for decades using quantitative data and
measurement tools [4, 5, 6]. Traditionally, meteorologists use various tools
for various purposes [7, 8, 9, 10]. Some of these tools include thermometers
for temperature measurement, barometers for air pressure measurement,
anemometers for wind speed measurement, etc. [11, 12, 13, 14]. With the
advancement of new tools and technologies and an increase in data and
computing power, we can use the data collected by these devices to predict the
weather using modern techniques which include statistical techniques and
machine learning [15, 16, 17, 18]. Weather prediction can be called a form of
data mining that is concerned with finding hidden patterns inside largely
available meteorological data [19, 20, 21, 22].

An accurate weather prediction is advantageous for a significant amount of the
population. In fact, every person and everything happening around us is
affected by the current weather condition [23, 24, 25, 26]. Some of the areas
include disaster management farming, navigation for ships and airplanes,
operation of hydro-power plants, sports and entertainment events, and even
common household activities [27, 28, 29, 30]. The process of trying to predict
the weather based on recurring meteorological and astronomical events began
with early human civilizations [31, 32, 33].

A vast amount of meteorological data is available even for the general public
which encourages the need to monitor the weather and atmospheric changes on a
daily basis [34, 35, 36, 37]. Many satellites just record atmospheric
parameters and conditions [38, 39, 40, 41]. Thus this large amount of data can
be used by machine learning models as an advantage over traditional techniques
to save time and increase accuracy [42, 43, 44, 45]. We used Multiple Linear
Regression with the data that are given at hand to predict the temperature of
a city for the upcoming days. The main focus is made on accuracy which has to
be more than the traditional methods.

##  2 Related Work

As the data is available cheaply and in vast quantity, researchers have
previously tried to propose models for predicting the weather. The nature of
weather is of course non-linear which makes the accuracy of prediction fairly
lower as expected [46, 47, 48, 49, 50]. In paper [51], a review of scientific
studies based on the prediction of weather using artificial neural networks is
carried out. The benefits of using artificial neural networks consisted of
yielding good results. Hence it can be used as an easy alternative to the
traditional meteorological approach. Artificial Neural Network can be seen as
capable of predicting most of the weather phenomena including rainfall, wind
speed, temperature, etc. [52] used backpropagation neural network model and
tested the idea with a real-time dataset. They compared the results with the
practical work of the meteorological department. The results were confirmed
that the model has the potential for successful application to temperature
forecasting.

Predict and classification of thunderstorms using ANN is also done by [53].
They designed the model to predict the occurrence of thunderstorms in two
geographical regions. ANN can be used effectively for the forecasting of
thunderstorms with a more than satisfactory level of accuracy. Derivation of
artificial neural networks for the purpose of weather prediction for a
particular location is done by [54]. They used a backpropagation neural
network model for initial modeling. Then the results obtained by the BPN model
are used to feed the Hopfield Networks. They used attributes such as wind
speed, humidity, and temperature. They collected three years of data
comprising of around 15000 instances.[55] applied the K- Nearest Neighbor
method to find the hidden pattern in a large dataset containing meteorological
data of a particular area. They achieved a very high rate of accuracy but
their model is unsuitable to use in remote areas.

##  3 Proposed Model

The working of the proposed model is described in this section followed by
accuracy computation of the model. We are trying to predict the mean
temperature of a day using the weather data for the last 3 days and the
regression technique of machine learning. We have used the Multiple Linear
Regression algorithm to train our model. Also, we have applied the feature
selection method to select the most decisive features out of all the features
available in our initial dataset. A diagram of this model can be seen in Fig.
1.

Figure 1: Block diagram of proposed model

The initial dataset is passed to the feature selection process which does two
tasks, first, it removes all the features that are not linearly dependent on
the dependent variable, and then it removes all the features that are not
significant. The output of this block is a refined dataset which is then
passed to the next process which uses this dataset to train our machine. The
output of the training process is a continuous function that is used to
predict the temperature.

###  3.1 Feature Selection

An important assumption in applying multiple linear regression is that there
should be a linear relationship between dependent and independent variables.
To check the linear behavior between dependent and independent variables, we
calculate the Pearson correlation coefficient ‘r’between all possible pairs of
dependent and independent variables. Pearson’s correlation coefficient’s value
lies between -1 and 1 inclusively, where negative value shows negative
correlation and positive value shows positive correlation. For a sufficient
linear relationship to exist, the value of r should be either less than -0.6
or more than 0.6. So, we discarded all the features where the absolute value
of r is less than 0.6.

After this, we applied the Backward Elimination method to remove insignificant
features. In Backward elimination, we use a hypothesis test. The idea is to
check the effect of the removal of a feature on the model. We calculate the
P-value to prove the insignificance of a feature. We first select a
significance level (in our case it is 5%), then we train our model and check
the P-value of all the features. The maximum P-value found is compared with
the significant level and if found greater, the corresponding feature is
discarded and the process is repeated with the remaining features until all
the features have a P-value lesser than the significant level.

###  3.2 Model Training by Multiple Linear Regression

After feature selection is done, we obtain a refined dataset. We use this
dataset to train our machine using a multiple linear regression algorithm. The
dataset is split into two parts, one is used to train the model and the other
is used to test the performance of the prediction. The training set stores 80%
of the data and the test set stores 20% of the data.

The formula to predict the value in multiple linear regression model (or the
regression function) is given in Eq. (1). Here, $k_{0}$, $k_{1}$….$k_{n}$ in
are the free parameters, y is the predicted value, and $x_{1}$, $x_{2}$….
$x_{n}$ are the values of the features. Once the value is predicted, we
calculate the square of the difference between actual and predicted values.
Then we add all the square values, and this sum is known as the error. The
objective of the algorithm is to minimize this error.

$y=k_{0}+{k_{1}}*x_{1}+{k_{2}}*x_{2}+{k_{3}}*x_{3}.......+{k_{n}}*x_{n}$ (1)

##  4 Results and Analysis

###  4.1 Experimental Setup

For our model, the dataset that we have used is created using the Weather
Underground’s API web service. This dataset has 997 instances, and each
instance stores the mean temperature of a day along with the weather details
for the last 3 days. The weather details for the last three days include mean
temperature, mean dew point temperature, mean humidity, precipitation, etc. of
each day.

| meantempm
---|---
maxdewptm_3 | 0.829230
maxtempm_3 | 0.832974
mindewptm_3 | 0.833546
meandewptm_3 | 0.834251
mimtempm_3 | 0.836340
maxdewptm_2 | 0.839893
meamdewptm_2 | 0.848907
mindewptm_2 | 0.852760
mintempm_2 | 0.854320
meamtempm_3 | 0.855662
maxtempm_2 | 0.863906
meantempm_2 | 0.881221
maxdewptm_1 | 0.887235
meandewptm_1 | 0.896681
mindewptm_1 | 0.899000
mintempm_1 | 0.905423
maxtempm_1 | 0.923787
meantempm_1 | 0.937563
meantempm | 1.000000
Figure 2: Absolute values of correlation coefficient of all the possible pairs
of dependent and independent variables sorted in ascending order

###  4.2 Feature Selection Results

To select linear features, the ‘corr’function of Pandas library in Python was
used to calculate the correlation coefficient(r). Fig. 2 shows the value of
$r$ after the removal of non-linear features. You can see in the diagram that
all the values are greater than 0.6. We have used the Stats Models Library
available in Python to get the P-value for all features. Fig. 3 shows the
results of the final features selection. Here, it can be seen that all the
P-values are less than the significant level (0.5). So, after the completion
of the feature selection step, we are left with 7 features that will be used
as input to multiple linear regression.

Figure 3: Summary of trained model after applying backward elimination

###  4.3 Performance Results

We used sklearn to train our model. Sklearn is a library for python users, it
provides multiple methods that implement various machine learning algorithms
and other different machine learning tasks. After training our model, we
predict the temperature of test set instances, we then plot a graph using
matplotlib in Python to have a graphical view of the accuracy of our model. We
draw a scatter plot between predicted values and actual values. The closer the
scatter plot is to line y = x, the model is more accurate. Fig. 4 shows the
scatter plot of the results obtained. We can see in Fig. 3 that our scatter
plot is closer to line y = x, which implies good accuracy.

Figure 4: Scattered plot diagram between actual and predicted values

##  5 Conclusion

A significant number of study has been done on weather forecasting in the past
to benefit human society as an accurate weather prediction can contribute to
saving losses including financial and human losses, and to predict day-to-day
functioning. Meteorological data is available freely or at a low cost with
numerous features and instances. Feature selection is important to reduce the
number of redundant and unrequired features. We started with more than 30
features and they were reduced to 7. We calculated the absolute mean error on
our test set. The absolute mean error thus obtained in the end is 2.8, which
implies that our model can predict the mean temperature of a day with an error
of 2.8 degrees Celsius, given the weather information of the last three days.
The deviation of 2.8 degrees Celsius is not much so the model can be used in a
simulation to get the idea of what the temperature is going to be.

## References

  * [1] D. Saxena, I. Gupta, J. Kumar, A. K. Singh, and X. Wen, “A Secure and Multiobjective Virtual Machine Placement Framework for Cloud Data Center,” _IEEE Systems Journal_ , pp. 1–12, 2021.
  * [2] I. Gupta and A. K. Singh, “GUIM-SMD: Guilty User Identification Model using Summation Matrix-based Distribution,” _IET Information Security_ , vol. 14, pp. 773–782, November 2020.
  * [3] A. Acharya, H. Prasad, V. Kumar, I. Gupta, and A. K. Singh, “MACI: Malicious API Call Identifier Model to Secure the Host Platform,” in _Proceedings of the Seventh International Conference on Mathematics and Computing_. Singapore: Springer Singapore, 2022, pp. 309–320.
  * [4] “Weather Forecasting – Wikipedia,” 2022. [Online]. Available: https://en.wikipedia.org/wiki/Weather_forecasting
  * [5] I. Gupta, R. Gupta, A. K. Singh, and R. Buyya, “MLPAM: A Machine Learning and Probabilistic Analysis Based Model for Preserving Security and Privacy in Cloud Environment,” _IEEE Systems Journal_ , vol. 15, no. 3, pp. 4248–4259, 2021.
  * [6] I. Gupta, “A Comparative Study of the Approach Provided for Preventing the Data Leakage,” _Other Topics Engineering Research eJournal_ , vol. 9, no. 5, September 2017.
  * [7] I. Gupta and A. K. Singh, “SELI: Statistical Evaluation based Leaker Identification Stochastic Scheme for Secure Data Sharing,” _IET Communications_ , vol. 14, pp. 3607–3618, December 2020.
  * [8] P. Godha, S. Jadon, A. Patle, I. Gupta, B. Sharma, and A. K. Singh, “Flooding and Forwarding Based on Efficient Routing Protocol,” in _International Conference on Innovative Computing and Communications_ , vol. 1166. Singapore: Springer Singapore, 2021, pp. 215–223, advances in Intelligent Systems and Computing.
  * [9] I. Gupta and A. K. Singh, “A Confidentiality Preserving Data Leaker Detection Model for Secure Sharing of Cloud Data using Integrated Techniques,” in _2019 7th International Conference on Smart Computing Communications (ICSCC)_. Curtin University, Sarawak Malaysia: IEEE, 2019, pp. 1–5.
  * [10] P. Tiwari, S. Mehta, N. Sakhuja, I. Gupta, and A. K. Singh, “Hybrid Method in Identifying the Fraud Detection in the Credit Card,” in _Evolutionary Computing and Mobile Sustainable Networks_ , vol. 53. Singapore: Springer Singapore, 2021, pp. 27–35, data Engineering and Communications Technologies.
  * [11] I. Gupta and A. K. Singh, “A Holistic View on Data Protection for Sharing, Communicating, and Computing Environments: Taxonomy and Future Directions,” 2022\.
  * [12] A. Kesharwani, A. Nag, A. Tiwari, I. Gupta, B. Sharma, and A. K. Singh, “Real-Time Human Locator and Advance Home Security Appliances,” in _Evolutionary Computing and Mobile Sustainable Networks_ , vol. 53. Singapore: Springer Singapore, 2021, pp. 37–49, data Engineering and Communications Technologies.
  * [13] K. Kaur, I. Gupta, and A. K. Singh, “A Comparative Study of the Approach Provided for Preventing the Data Leakage,” _International Journal of Computer Networks & Communications (IJCNC)_.
  * [14] G. Batra, H. Singh, I. Gupta, and A. K. Singh, “Best Fit Sharing and Power Aware (BFSPA) Algorithm for VM Placement in Cloud Environment,” in _2017 3rd International Conference on Advances in Computing,Communication & Automation (ICACCA) (Fall)_. IEEE, 2017, pp. 1–4.
  * [15] “Students Ask: How Do Meteorologists Predict The Weather,” February 2012. [Online]. Available: https://www.gpb.org/blogs/talking-up-a-storm/2012/02/15/students-ask-how-do-meteorologists-predict-the-weather
  * [16] I. Gupta, N. Singh, and A. Singh, “Layer-based Privacy and Security Architecture for Cloud Data Sharing,” _Journal of Communications Software and Systems (JCOMSS)_ , vol. 15, no. 2, 2019.
  * [17] I. Gupta and K. Gupta, “Review on Intrusion Detection System Architectures in WSN,” _International Journal of Scientific & Engineering Research_, vol. 7, no. 12, pp. 111–115, Dec. 2016.
  * [18] K. Kaur, I. Gupta, and A. K. Singh, “A Comparative Evaluation of Data Leakage/Loss Prevention Systems (DLPS),” in _Proc. 4th International Conference Computer Science & Information Technology_, 2017, pp. 87–95.
  * [19] S. N. Kohail and A. M. El-Halees, “Implementation of Data Mining Techniques for Meteorological Data Analysis,” _International Journal of Information and Communication Technology Research_ , vol. 1, no. 3, pp. 96–100, 2011.
  * [20] V. Sharma, S. Jalwa, A. R. Siddiqi, I. Gupta, and A. K. Singh, “A Lightweight Effective Randomized Caesar Cipher Algorithm for Security of Data,” in _Evolutionary Computing and Mobile Sustainable Networks_ , vol. 53. Singapore: Springer Singapore, 2021, pp. 411–419, data Engineering and Communications Technologies.
  * [21] U. Arora, S. Verma, I. Gupta, and A. K. Singh, “Implementing Privacy using Modified Tree and Map Technique,” in _2017 3rd International Conference on Advances in Computing,Communication & Automation (ICACCA) (Fall)_. IEEE, 2017, pp. 1–5.
  * [22] Khushbu, P. Nishad, V. Kashyap, and I. Gupta, “A Classification and Distribution Model for Data Leakage Prevention and Detection,” _International Research Journal of Modernization in Engineering Technology and Science_ , vol. 3, no. 2, pp. 348–354, Feb. 2021.
  * [23] I. Gupta and A. K. Singh, “Dynamic Threshold based Information Leaker Identification Scheme,” _Information Processing Letters_ , vol. 147, pp. 69 – 73, 2019.
  * [24] A. Acharya, H. Prasad, V. Kumar, I. Gupta, and A. K. Singh, “Host Platform Security and Mobile Agent Classification: A Systematic Study,” in _Computer Networks and Inventive Communication Technologies_ , vol. 58. Singapore: Springer Singapore, 2021, pp. 1001–1010, data Engineering and Communications Technologies.
  * [25] P. Agarwal, S. Mittal, A. Tiwari, I. Gupta, A. K. Singh, and B. Sharma, “Authenticating Cryptography over Network in Data,” in _2019 International Conference on Intelligent Computing and Control Systems (ICCS)_. IEEE, 2019, pp. 632–636.
  * [26] A. K. Singh, I. Gupta, R. Verma, V. Gautam, and C. P. Yadav, “A Survey on Data Leakage Detection and Prevention,” in _Proc. Int. Conf. Innov. Comput. Commun._ , 2020.
  * [27] I. Gupta and A. K. Singh, “An Integrated Approach for Data Leaker Detection in Cloud Environment,” _Journal of Information Science and Engineering_ , vol. 36, pp. 993–1005, Sep. 2020.
  * [28] R. Verma, V. Gautam, C. P. Yadav, I. Gupta, and A. K. Singh, “A Survey on Data Leakage Detection and Prevention,” in _Proceedings of the International Conference on Innovative Computing & Communications (ICICC) 2020_. SSRN, Elsevier, May 2020, pp. 1–7.
  * [29] S. Jalwa, V. Sharma, A. R. Siddiqi, I. Gupta, and A. K. Singh, “Comprehensive and Comparative Analysis of Different Files Using CP-ABE,” in _Advances in Communication and Computational Technology_ , vol. 668. Singapore: Springer Singapore, 2021, pp. 189–198, electrical Engineering.
  * [30] Khushbu, P. Nishad, V. Kashyap, I. Gupta, and A. K. Singh, “An Organized Study on Data Divulge Elimination and Discernment,” in _Computer Networks and Inventive Communication Technologies_. Singapore: Springer Singapore, 2021, pp. 569–578.
  * [31] D. Chauhan and J. Thakur, “Data Mining Techniques for Weather Prediction: A Review,” _International Journal on Recent and Innovation Trends in Computing and Communication_ , vol. 2, no. 8, pp. 2184–2189, 2014.
  * [32] I. Gupta and K. Gupta, “Evaluation of Intrusion Detection Schemes in Wireless Sensor Network,” _IOSR Journal of Computer Engineering_ , vol. 18, no. 2, pp. 60–63, Mar-Apr. 2016.
  * [33] I. Gupta and A. K. Singh, “A Probabilistic Approach for Guilty Agent Detection using Bigraph after Distribution of Sample Data,” _Procedia Computer Science_ , vol. 125, pp. 662 – 668, 2018.
  * [34] K. Kaur, I. Gupta, and A. K. Singh, “A Comparative Study of the Approach Provided for Preventing the Data Leakage,” _International Journal of Network Security & Its Applications_, vol. 9, no. 5, pp. 21–33, 2017.
  * [35] A. Nag, A. Kesharwani, B. Sharma, I. Gupta, A. Tiwari, and A. K. Singh, “Potential and Extention of Internet of Things,” in _Second International Conference on Computer Networks and Communication Technologies (ICCNCT)_ , vol. 44. Cham: Springer International Publishing, 2020, pp. 542–551.
  * [36] K. Kaur, I. Gupta, and A. K. Singh, “E-Mail Protection System to Prevent Data Leakage,” _Vigyan Prakash_ , vol. 16, pp. 30–36, 2018.
  * [37] K. N. Kaur, Divya, I. Gupta, and A. K. Singh, “Digital Image Watermarking Using (2, 2) Visual Cryptography with DWT-SVD Based Watermarking,” in _Computational Intelligence in Data Mining_ , vol. 711. Singapore: Springer Singapore, 2019, pp. 77–86, advances in Intelligent Systems and Computing.
  * [38] I. Gupta and A. K. Singh, “A Framework for Malicious Agent Detection in Cloud Computing Environment,” _International Journal of Advanced Science and Technology (IJAST)_ , vol. 135, pp. 49–62, Feb 2020.
  * [39] K. Kaur, I. Gupta, and A. K. Singh, “Data Leakage Prevention: E-Mail Protection via Gateway,” _Journal of Physics: Conference Series_ , vol. 933, p. 012013, jan 2018.
  * [40] K. Gupta and I. Gupta, “A Comprehensive Study on Architecture, Security issues and Challenges in Cloud Computing,” _International Journal of Scientific & Engineering Research_, vol. 7, no. 12, pp. 128–131, Dec. 2016.
  * [41] P. Godha, S. Jadon, A. Patle, I. Gupta, B. Sharma, and A. Kumar Singh, “Architecture, an Efficient Routing, Applications, and Challenges in Delay Tolerant Network,” in _2019 International Conference on Intelligent Computing and Control Systems (ICCS)_. IEEE, 2019, pp. 824–829.
  * [42] E. Abrahamsen, O. Brastein, and B. Lie, “Machine Learning in Python for Weather Forecast based on Freely Available Weather Data,” in _59th Conference on Simulation and Modelling (SIMS 59)_ , 11 2018, pp. 169–176.
  * [43] V. Gautam, C. P. Yadav, R. Verma, I. Gupta, and A. K. Singh, “A Comparative Study on Data Leakage Detection and Prevention,” in _IEEE International Conference on Nanoelectronics, Nanophotonics, Nanomaterials, Nanobioscience & Nanotechnology (5NANO2021)_. IEEE, 2021.
  * [44] I. Gupta and A. K. Singh, “A Hybrid Technique for the Detection of Data Leakage in Cloud computing Environment,” in _Ist International Conference on Science in Hindi_ , August 2017, vigyan Prakash.
  * [45] A. K. Singh and I. Gupta, “Online Information Leaker Identification Scheme for Secure Data Sharing,” _Multimedia Tools and Applications_ , vol. 79, no. 41, pp. 31 165–31 182, November 2020.
  * [46] I. Gupta and A. K. Singh, “A Probability based Model for Data Leakage Detection using Bigraph,” in _Proceedings of 7th International Conference on Communication and Network Security (ICCNS)_ , ser. ICCNS 2017\. New York, NY, USA: Association for Computing Machinery (ACM), 2017, p. 1–5.
  * [47] I. Gupta, T. K. Madan, S. Singh, and A. K. Singh, “HISA-SMFM: Historical and Sentiment Analysis Based Stock Market Forecasting Model,” 2022.
  * [48] I. Gupta, V. Sharma, S. Kaur, and A. K. Singh1, “PCA-RF: An Efficient Parkinson’s Disease Prediction Model based on Random Forest Classification,” 2022.
  * [49] I. Gupta, S. Mittal, A. Tiwari, P. Agarwal, and A. K. Singh, “TIDF-DLPM: Term and Inverse Document Frequency based Data Leakage Prevention Model,” 2022.
  * [50] I. Gupta, P. K. Yadav, S. Pareek, S. Shakeel, and A. K. Singh, “Auxiliary Informatics System: an Advancement towards a Smart Home Environment,” 2022.
  * [51] A. Saxena, N. Verma, and K. Tripathi, “A Review Study of Weather Forecasting Using Artificial Neural Network Approach,” _International Journal of Engineering Research & Technology_, vol. 2, no. 11, pp. 2029–2035, 2013.
  * [52] J. Devi, B. Reddy, K. Kumar, B. Reddy, and N. RajaNayak, “ANN Approach for Weather Prediction Using Back Propagation,” _International Journal of Engineering Trends and Technology_ , vol. 3, no. 1, pp. 19–23, 2012.
  * [53] S. M. Anad, A. Dash, M. S. J. Kumar, and A. Kesarkar, “Prediction and Classification of Thunderstorms Using Artificial Neural Network,” _International Journal of Engineering Science and Technology_ , vol. 3, no. 5, pp. 4031–4035, 2011.
  * [54] G. J. Sawale and S. R. Gupta, “Use of Artificial Neural Network in Data Mining For Weather Forecasting,” _International Journal Of Computer Science And Applications_ , vol. 6, no. 2, pp. 384–387, 2013.
  * [55] S. Badhiye, W. V, and C. N, “Analysis of Temperature and Humidity Data for Future Value Prediction,” _International Journal of Computer Science and Information Technologies_ , vol. 3, no. 1, pp. 3012–3014, 05 2012.