# Multiscale Modeling and Analysis of Growth of Plant Tissues

Arezki Boudaoud, Annamaria Kiss, Mariya Ptashnyk

► To cite this version:

Arezki Boudaoud, Annamaria Kiss, Mariya Ptashnyk. Multiscale Modeling and Analysis of Growth of Plant Tissues. SIAM Journal on Applied Mathematics, 2023, 83 (6), pp.2354-2389. 10.1137/23M1553315. hal-04334494

HAL Id: hal-04334494

https://hal.science/hal-04334494v1

Submitted on 30 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
MULTISCALE MODELLING AND ANALYSIS OF GROWTH OF PLANT TISSUES

AREZKI BOUDAOUD, ANNAMARIA KISS, MARIYA PTASHNYK

LADHYX, CNRS, ECOLE POLYTECHNIQUE, IP PARIS, PALAISEAU, FRANCE

RDP, ENS DE LYON, CLAUDE BERNARD UNIVERSITY LYON 1, CNRS, INRAE, LYON, FRANCE

DEPARTMENT OF MATHEMATICS, HERIOT-WATT UNIVERSITY, THE MAXWELL INSTITUTE FOR MATHEMATICAL SCIENCES, EDINBURGH, SCOTLAND, UK

**ABSTRACT.** How morphogenesis depends on cell properties is an active direction of research. Here, we focus on mechanical models of growing plant tissues, where microscopic (sub)cellular structure is taken into account. In order to establish links between microscopic and macroscopic tissue properties, we perform a multiscale analysis of a model of growing plant tissue with subcellular resolution. We use homogenization to rigorously derive the corresponding macroscopic tissue scale model. Tissue scale mechanical properties are computed from microscopic structural and material properties, taking into account deformation by the growth field. We then consider case studies and numerically compare the detailed microscopic model and the tissue-scale model, both implemented using finite element method. We find that the macroscopic model can be used to efficiently make predictions about several configurations of interest. Our work will help making links between microscopic measurements and macroscopic observations in growing tissues.

# 1. INTRODUCTION

Modelling plant growth and morphogenesis is an active area of research [53]. A major difficulty in this area is that growth is ‘inherently a multiscale process’ [53], bridging subcellular to organ scales. Each plant cell is surrounded by a thin layer of polysaccharides, known as the cell wall, and exerts on this wall a hydrodynamic pressure, termed turgor pressure. Plant cell growth is driven by turgor pressure and restrained by the cell wall. Organ morphogenesis results from specific spatial distributions of growth rates across cells and tissues. As a consequence, the form of an organ is dependent on all processes occurring from sub-wall scale to supra-cellular scale, raising the need for multiscale modelling approaches.

Most previous modelling effort has focused on one scale or one process at a time. For instance, several studies have addressed how cell wall rheology emerges from its composition, microstructure, and/or synthesis [7, 19, 20, 28, 44]. These studies adopted approaches from continuum mechanics and used partial differential equations (PDEs) to describe the cell wall as an anisotropic material, with different assumptions on the rheology – viscous [19, 20], elastic [7, 8, 44, 47], or viscoelastoplastic [28]; growth then corresponds to flow of the viscous material or to remodelling of the elastic material, combined with synthesis of new material. Other studies represented tissues as tessellations by polygons of 2D space or of surfaces embedded in 3D that describe the position of vertices, relying on high-dimensional systems of ordinary differential equations [7, 15, 24, 31, 35]. The authors investigated how tissue shape changes according to cell wall rheology, assuming for instance each edge to be a viscoelastoplastic element, or to how mechanical stress in tissues feeds back on cell wall dynamics. Other studies used continuous, PDE-based approaches for modelling tissue dynamics and investigated how patterns of tissue mechanical properties or of growth potential yield organ morphogenesis [26, 54], see the comprehensive overview in [22].

In this context, a systematic derivation of macroscopic (supra-cellular) tissue rheology from (sub)wall rheology is still lacking. Here we try to address this issue using multiscale modelling and homogenization techniques. There are many results on homogenization of equations of linear elasticity, e.g. [4, 14, 29, 43, 45, 50], however to our knowledge the multiscale analysis of the two-way coupling between elastic deformation and growth presented in here is novel.

In the derivation of the microscopic model we start from the framework of morphoelasticity [23, 49]. Following [49], we consider the multiplicative decomposition of the deformation gradient into elastic and growth part and model elastic deformations of plant cell walls using equations of linearised elasticity. To describe the growth we use Lockhart’s law [32] that relates plant cell wall growth with deformation gradient and accounts for microscopic subwall properties. Along with modelling growth processes, the multiplicative decomposition approach is also used in multiscale modelling and analysis of plastic deformations, where the decomposition corresponds to the elastic and plastic deformations respectively, e.g. [16, 36, 39]. Rigorous well-posedness results were recently obtained for a model combining multiplicative decomposition and nonlinear elasticity [17].
Applying homogenization techniques, the two-scale convergence [3, 41] and periodic unfolding [11, 13] methods, we rigorously derive macroscopic equations that describe elastic deformations and growth at the tissue scale. In particular, we provide explicit formulae for macroscopic elastic properties as a function of microscopic parameters and microscopic structure. An important step in the rigorous derivation of the macroscopic problem is the proof of the strong convergence of the sequence of growth and strain tensors as the small parameter, representing a ratio between the size of the microstructure and the size of the tissue, converges to zero. We illustrate our results by solving numerically the corresponding microscopic and macroscopic equations and analyse the degree of agreement between solutions of the macroscopic problem and solutions of the original microscopic model defined at the cell wall scale. An important contribution to the numerical simulation of the macroscopic two-scale problem is the development of a two-scale numerical algorithm that allows an efficient coupling between macroscopic and microscopic properties and processes. Along with numerical efficiency, the advantage of the derivation of macroscopic tissue level models for plant growth allow us to study different biological settings which are difficult or even impossible to formulate and simulate at the cell-scale level.

The paper is organized as follows. In section 2 we formulate and analyse the microscopic model for plant tissue growth. The macroscopic model is derived in sections 3 using both formal asymptotic expansion and rigorous derivation applying the two-scale convergence and the periodic unfolding methods. Description of the numerical simulation algorithm and implementation of the numerical methods is given in section 4. In section 5 we present numerical simulation results, followed by a short conclusion and an appendix.

## 2. DERIVATION OF MICROSCOPIC MODEL FOR PLANT TISSUE GROWTH

In our mathematical model for the growth of a plant tissue we consider a microscopic geometry of a plant tissue composed of cells surrounded by cell walls, connected by middle lamella, see specific geometry in Figure 1a-b. In modelling the dynamics of plant tissue we consider the elastic deformations and growth of these cell walls and middle lamella.

By $\Omega_t \subset \mathbb{R}^d$, with $t > 0$ and $d = 2$ or 3, we shall represent a part of a plant tissue in the current configuration at time $t$ and $\partial\Omega_t$ denotes the external boundary of the tissue. We shall consider the growth and elastic deformation of a plant tissue given by the map $\chi(t, \cdot) : \Omega \to \mathbb{R}^d$ from the initial (reference) configuration $\Omega \subset \mathbb{R}^d$ into deformed (current) configuration $\Omega_t = \chi(t, \Omega)$ of a plant tissue for $t > 0$. We consider $\Omega$ to be a bounded Lipschitz convex or a bounded $C^{1,\tilde{\gamma}}$, with $\tilde{\gamma} \in (0, 1)$, domain.

When modelling growth we use the framework of morphoelasticity and consider the multiplicative decomposition of the deformation gradient $\mathbf{F} = \mathbf{I} + \nabla\mathbf{u}$ into elastic and growth parts $\mathbf{F} = \mathbf{F}_e\mathbf{F}_g$, where $\mathbf{u}$ is the displacement of the cell walls and middle lamella according to the map $\chi$ and $\mathbf{F}_e$ and $\mathbf{F}_g$ are elastic and growth deformation gradients, respectively. Then the elastic strain is given by

$$ (1) \qquad \begin{aligned} \mathbf{E}^{\text{el}} &= \frac{1}{2}(\mathbf{F}_e^T \mathbf{F}_e - \mathbf{I}) = \frac{1}{2}[(\mathbf{F}_e - \mathbf{I})^T(\mathbf{F}_e - \mathbf{I}) + (\mathbf{F}_e - \mathbf{I}) + (\mathbf{F}_e - \mathbf{I})^T] \\ &= \frac{1}{2}[(\mathbf{F}_e - \mathbf{I})^T(\mathbf{F}_e - \mathbf{I}) + (\mathbf{I} + \nabla\mathbf{u})\mathbf{F}_g^{-1} + ((\mathbf{I} + \nabla\mathbf{u})\mathbf{F}_g^{-1})^T - 2\mathbf{I}]. \end{aligned} $$

We model the cell walls and middle lamella as an hyperelastic material and consider the stress tensor in the form

$$ \sigma(x, \mathbf{F}_e) = J_e^{-1} \mathbf{F}_e \partial_{\mathbf{F}_e} W(x, \mathbf{F}_e)^T, $$

where $W$ is the strain energy function and $J_e = \det(\mathbf{F}_e) = \det(\mathbf{F}\mathbf{F}_g^{-1})$. Then the constitutive equation for elastic deformations is given by

$$ \sigma(x, \mathbf{F}_e) = \sigma(x, \nabla\mathbf{u}, \mathbf{F}_g^{-1}) = J_e^{-1} \mathbf{F}\mathbf{F}_g^{-1} \partial_{\mathbf{F}_e} W(x, \mathbf{F}\mathbf{F}_g^{-1})^T. $$

Mechanical equilibrium requires that

$$ (2) \qquad -\text{div}_{\tilde{x}} \sigma(x, \nabla\mathbf{u}, \mathbf{F}_g^{-1}) = 0 \quad \text{in } \Omega_t^w, \ t > 0, $$

where $\Omega_t^w$ denotes the domain of cell walls and middle lamellae, joining the walls of individual cells together, and $\tilde{x}$ denotes the coordinates in the current configuration. We complete equations (2) with the boundary conditions

$$ (3) \qquad \begin{aligned} \sigma(x, \nabla\mathbf{u}, \mathbf{F}_g^{-1})\nu &= -P\nu && \text{on } \Gamma_t, \ t > 0, \\ \sigma(x, \nabla\mathbf{u}, \mathbf{F}_g^{-1})\nu &= \mathbf{f} && \text{on } \partial\Omega_t \setminus \Gamma_{D,t}, \ t > 0, \\ \mathbf{u} \cdot \nu &= 0, \quad \Pi_{\tau}(\sigma(x, \nabla\mathbf{u}, \mathbf{F}_g^{-1})\nu) = 0 && \text{or } \mathbf{u} = 0 && \text{on } \Gamma_{D,t}, \ t > 0, \end{aligned} $$

where $P$ is the turgor pressure inside the cells which can vary between cells and across the tissue, $\mathbf{f}$ denotes external forces, $\nu$ is the external normal vector to the boundaries of $\Omega_t^w$ in current configuration, $\Gamma_t$ denotes the
boundary of cells, corresponding to plasma membranes, $\Gamma_{D,t} \subset \partial\Omega_t$ for $t \ge 0$, and $\Pi_{\tau}$ denotes the tangential components of the corresponding vector.

We can rewrite (2) and (3), defined in the current configuration, in the reference configuration to obtain

$$ (4) \qquad \begin{aligned} -\operatorname{div}(J_g \mathbf{S}^T(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \mathbf{F}_g^{-T}) &= 0 && \text{in } \Omega^w, t > 0, \\ J_g \mathbf{S}^T(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \mathbf{F}_g^{-T} N &= -J P \mathbf{F}^{-T} N && \text{on } \Gamma, t > 0, \\ J_g \mathbf{S}^T(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \mathbf{F}_g^{-T} N &= \mathbf{f} J |\mathbf{F}^{-T} N| && \text{on } \partial\Omega \setminus \Gamma_D, t > 0, \\ \mathbf{u} \cdot N = 0, \quad \Pi_{\tau}(J_g \mathbf{S}^T(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \mathbf{F}_g^{-T} N) &= 0 \quad \text{or} \quad \mathbf{u} = 0 && \text{on } \Gamma_D, t > 0, \end{aligned} $$

where $N$ denotes the external normal vector to the boundaries of the reference domain of cell walls and middle lamella $\Omega^w \subset \mathbb{R}^d$, and the nominal elastic stress is given by

$$ \mathbf{S}^T(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) = \partial_{\mathbf{F}_e} W(x, \mathbf{F} \mathbf{F}_g^{-1}). $$

To specify the constitutive relation for the stress tensor $\sigma$, or correspondent nominal tensor $S$, we assume small elastic strain (small elastic deformations of plant tissues), i.e.

$$ \sigma(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \approx \mathbf{S}^T(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \approx \mathbb{E}(x) \varepsilon^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g^{-1}), $$

where $\mathbb{E}(x) = \partial_{\mathbf{F}_e} \partial_{\mathbf{F}_e} W(x, \mathbf{1})$ is the elasticity tensor and $\varepsilon^{\text{el}}$ is the linearised version of the elastic strain (1), which depends on the displacement gradient and growth tensor

$$ \varepsilon^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g^{-1}) = \frac{1}{2} [\nabla \mathbf{u} \mathbf{F}_g^{-1} + (\nabla \mathbf{u} \mathbf{F}_g^{-1})^T + \mathbf{F}_g^{-1} + \mathbf{F}_g^{-T} - 2\mathbf{I}] = \operatorname{sym}(\nabla \mathbf{u} \mathbf{F}_g^{-1}) + \operatorname{sym}(\mathbf{F}_g^{-1}) - \mathbf{I}. $$

We also assume $J\mathbf{F}^{-T} \approx J_g\mathbf{F}_g^{-T}$ for small elastic strain. Notice that for $\mathbf{F}_g = \mathbf{I}$ we recover the standard formula for the strain in the case of linear elasticity.

To complete the model we specify equation for the growth tensor $\mathbf{F}_g$

$$ (5) \qquad \begin{aligned} \frac{\partial \mathbf{F}_g}{\partial t} &= G(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \mathbf{F}_g && \text{in } \Omega^w, t > 0, \\ \mathbf{F}_g(0, x) &= \mathbf{I} && \text{for } x \in \Omega^w. \end{aligned} $$

Since there is no consensus on modelling growth [22, 53], we consider two scenarios and assume that the growth depends on the local average of stress or of strain in cell walls and middle lamella (both are compatible with Lockhart's law [32]). We also assume that the cell wall and middle lamella expand when the local average of the stress or strain is larger than some threshold value [32]. Hence we consider the stress based growth

$$ \tilde{G}(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) = \tilde{G}(x, \sigma(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1})) = \eta_{\sigma}[\hat{\sigma} - \tau_{\sigma}]_{+}, $$

or the elastic strain based growth

$$ \tilde{G}(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) = \tilde{G}(x, \varepsilon^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g^{-1})) = \eta_{\varepsilon}[\hat{\varepsilon}^{\text{el}} - \tau_{\varepsilon}]_{+}, $$

and

$$ (6) \qquad G(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1})_{ij} = \begin{cases} -M & \text{if } \tilde{G}(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1})_{ij} \le -M, \\ \tilde{G}(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1})_{ij} & \text{if } -M < \tilde{G}(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1})_{ij} < M, \\ M & \text{if } \tilde{G}(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1})_{ij} \ge M, \end{cases} \quad \text{for } i, j = 1, \dots, d, $$

for some $M > 0$. Here $[v]_+ = Q^T[Dv]^+Q$ with $Q$ the rotation that diagonalizes the tensor $v$, and $[Dv]^+$ is a diagonal matrix with the positive parts of the eigenvalues of $v$ on its diagonal. This allows to relate growth only to tensile stress, in the case of stress based growth, or only to elongational strain, in the case of strain based growth. The uniform boundedness assumption on the growth rate $G$ is used in the rigorous analysis of the model and is not restrictive from the biological point of view. Notice that in the growth laws we consider piece-wise constant fields $\hat{\sigma}$ and $\hat{\varepsilon}^{\text{el}}$, obtained as an average over each cell of $\sigma$ and $\varepsilon^{\text{el}}$ respectively. The growth laws introduce two physical parameters: the extensibility constants $\eta_{\sigma}$ and $\eta_{\varepsilon}$ and threshold matrices $\tau_{\sigma}$ and $\tau_{\varepsilon}$, respectively.
2.1. **Formulation of the microscopic model.** We assume that in a plant tissue cells are distributed periodically and consider the parameter $\delta > 0$ that determines the ratio between the size of a cell and the size of the tissue. We also assume that the size of the cell and the thickness of the cell wall $h$ are of the same order and much smaller than the size of the tissue, i.e. $\delta$ is small. To define the microscopic structure of the plant tissue given by the cell walls and middle lamella, we consider a 'unit cell' $Y$ and $Y_c \subset Y$ represents the voids filled by biological cells, with Lipschitz boundary $\Gamma = \partial Y_c$ and composed of a finite number of subdomains separated from each other and from the edges of $Y$, whereas $Y_w = Y \setminus \bar{Y}_c$ represents cell walls surrounded by middle lamella. Then the microscopic geometry of a plant tissue in the reference configuration is defined as

$$ \Omega_c^\delta = \bigcup_{\xi \in \Xi_\delta} \delta(\bar{Y}_c + \xi) \quad \text{and} \quad \Omega^\delta = \Omega \setminus \Omega_c^\delta, $$

where $\Xi_\delta = \{\xi \in \Xi : \delta(\bar{Y} + \xi) \subset \bar{\Omega}, \text{dist}(\delta(\bar{Y}_c + \xi), \partial\Omega) \ge \kappa\delta\}$ and $\Xi = \{\xi \in \mathbb{R}^d : \xi = \sum_{j=1}^d k_j b_j, k \in \mathbb{Z}^d\}$, with $\{b_j\}_{j=1}^d$ being the basis vectors of $Y$, i.e. $Y = \{y \in \mathbb{R}^d : y = \sum_{j=1}^d s_j b_j, s \in (0,1)^d\}$, and for some fixed $\kappa > 0$. The boundaries of $\Omega^\delta$ that correspond to cell plasma-membranes are

$$ \Gamma^\delta = \bigcup_{\xi \in \Xi_\delta} \delta(\Gamma + \xi) \quad \text{and} \quad \Gamma^\delta = \partial\Omega^\delta \setminus \partial\Omega. $$

The parts of the tissues near the boundary $\partial\Omega$, which do not include the complete $\delta\bar{Y}$ are denoted by

$$ \Lambda_{\delta} = \Omega \setminus \bigcup_{\xi \in \Xi_{\delta}} \delta(\bar{Y} + \xi) \quad \text{and} \quad \Lambda_{\delta} = \bigcup_{l=1}^{L} \Lambda_{\delta}^{l}, \quad \text{with} \quad \Lambda_{\delta}^{l} \cap \delta(Y + \xi) \neq \emptyset \quad \text{for one } \xi \in \Xi \text{ and } L = \mathcal{O}(1/\delta^{d-1}). $$

Then microscopic equations for elastic deformations in the reference domain read

$$ \begin{array}{ll} -\text{div}(J_g^\delta \mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) \mathbf{F}_{g,\delta}^{-T}) = 0 & \text{in } \Omega^\delta, t > 0, \\ J_g^\delta \mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) \mathbf{F}_{g,\delta}^{-T} N = -J_g^\delta P^\delta(t, x) \mathbf{F}_{g,\delta}^{-T} N & \text{on } \Gamma^\delta, t > 0, \\ J_g^\delta \mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) \mathbf{F}_{g,\delta}^{-T} N = J_g^\delta \mathbf{f}(t, x) |\mathbf{F}_{g,\delta}^{-T} N| & \text{on } \partial\Omega \setminus \Gamma_D, t > 0, \\ \mathbf{u}^\delta \cdot N = 0, \quad \Pi_\tau(J_g^\delta \mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) \mathbf{F}_{g,\delta}^{-T} N) = 0 \quad \text{or} \quad \mathbf{u}^\delta = 0 & \text{on } \Gamma_D, t > 0, \end{array} $$

where $J_g^\delta = \det(\mathbf{F}_{g,\delta})$, $\mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) = \mathbb{E}^\delta(x)\boldsymbol{\varepsilon}^{\text{el}}(\nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) = \mathbb{E}^\delta(x)[\text{sym}(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1}) + \text{sym}(\mathbf{F}_{g,\delta}^{-1}) - \mathbf{I}]$, with $\mathbb{E}^\delta(x) = \mathbb{E}(x, x/\delta)$, and $P^\delta(t, x) = P(t, x, x/\delta)$ for given functions $\mathbb{E}: \Omega \times Y \to \mathbb{R}^{4d}$ and $P: (0, T) \times \Omega \times \Gamma \to \mathbb{R}$, extended $Y$-periodically to $\mathbb{R}^d$ and to $\bigcup_{\xi \in \Xi}(\Gamma + \xi)$ respectively.

The dynamics of the growth tensor is determined by

$$ \begin{array}{ll} \displaystyle \frac{\partial \mathbf{F}_{g,\delta}}{\partial t} = G^\delta(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) \mathbf{F}_{g,\delta} & \text{in } \Omega^\delta, t > 0, \\ \mathbf{F}_{g,\delta}(0, x) = \mathbf{I} & \text{in } \Omega^\delta, \end{array} $$

where

$$ \begin{array}{l} \tilde{G}^\delta(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) = \eta_\sigma \left[ \int_{\delta([x/\delta]_Y + Y_w) \cap \Omega} \boldsymbol{\sigma}(\tilde{x}, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) d\tilde{x} - \boldsymbol{\tau}_\sigma \right]_+ \quad \text{or} \\ \tilde{G}^\delta(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) = \eta_\varepsilon \left[ \int_{\delta([x/\delta]_Y + Y_w) \cap \Omega} \boldsymbol{\varepsilon}^{\text{el}}(\nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) d\tilde{x} - \boldsymbol{\tau}_\varepsilon \right]_+ \end{array} $$

with $\boldsymbol{\sigma}(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) = \mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) = \mathbb{E}^\delta(x)\boldsymbol{\varepsilon}^{\text{el}}(\nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1})$, and

$$ (9) \quad G^\delta(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1})_{ij} = \begin{cases} -M & \text{if } \tilde{G}^\delta(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1})_{ij} \le -M, \\ \tilde{G}^\delta(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1})_{ij} & \text{if } -M < \tilde{G}^\delta(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1})_{ij} < M, \\ M & \text{if } \tilde{G}^\delta(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1})_{ij} \ge M, \end{cases} \quad \text{for } i, j = 1, \dots, d. $$

2.2. **Well-posedness of microscopic model.** To prove existence of a weak solution of problem (7) and (8), we consider standard ellipticity assumptions on the elasticity tensor $\mathbb{E}$ and regularity assumptions on the pressure $P^\delta$ and boundary forces $\mathbf{f}$. For the pressure inside of the cells $P^\delta$ we shall consider dependence on the microscopic structure in the form

$$ (10) \qquad P^\delta(t, x) = P_1(t, x) + \delta P_2(t, x, x/\delta), $$

for some given functions $P_1: (0, T) \times \Omega \to \mathbb{R}$ and $P_2: (0, T) \times \Omega \times \Gamma \to \mathbb{R}$.
**Assumption 2.1.** • Elasticity tensor $\mathbb{E} \in C^\gamma(\bar{\Omega}; L_\text{per}^\infty(Y))^{4d}$ is positive definite and bounded, i.e. $\alpha_1|\mathbf{A}|^2 \le \mathbb{E}(x, y)\mathbf{A} \cdot \mathbf{A} \le \alpha_2|\mathbf{A}|^2$ for $x \in \Omega$, $y \in Y$, symmetric matrices $\mathbf{A} \in \mathbb{R}^{d \times d}$, and positive constants $\alpha_1, \alpha_2$, and has minor and major symmetries, i.e. $\mathbb{E}_{ijkl} = \mathbb{E}_{klij} = \mathbb{E}_{jikl} = \mathbb{E}_{ijkl}$, for $i, j, k, l = 1, \dots, d$, for $\gamma \in (0, 1)$.

• $\mathbf{f} \in C^\gamma([0, T] \times \partial\Omega)$, $P_1 \in C^\gamma([0, T]; C^1(\bar{\Omega}))$ and $P_2 \in C^\gamma([0, T] \times \bar{\Omega} \times \Gamma)$, where $P_2(t, x, \cdot)$ is $Y$-periodically extended to $\Gamma + \xi$, with $\xi \in \Xi$, for $\gamma \in (0, 1)$.

If $\Gamma_D \neq \emptyset$, consider $\mathcal{G}$ the symmetry group of $\Gamma_D$, formed of at most one, for $d=2$, or two, for $d=3$, translations and at most one rotation, for $d=3$, that leave $\Gamma_D$ invariant. Let $\mathcal{G}_\tau$ be the subspace of $\mathbb{R}^d$ spanned by the set of translations in $\mathcal{G}$ and $\{\rho_1, \rho_2\}$ be the orthonormal basis of the plane perpendicular to the rotation axis. Then, depending on the boundary conditions, we define the following space for solutions of (7)

$$V_{\delta} = \{\mathbf{u} \in H^1(\Omega^{\delta})^d : \mathbf{u} = 0 \text{ on } \Gamma_D\} \quad \text{or} \\ V_{\delta} = \{\mathbf{u} \in H^1(\Omega^{\delta})^d : \mathbf{u} \cdot \mathbf{N} = 0 \text{ on } \Gamma_D, \int_{\Omega^{\delta}} \Pi_{\mathcal{G}_{\tau}}(\mathbf{u}) \, dx = 0, \int_{\Omega} \rho_1(\nabla \hat{\mathbf{u}} - \nabla \hat{\mathbf{u}}^T) \rho_2 \, dx = 0\},$$

where $\Pi_{\mathcal{G}_\tau}$ is the projection on $\mathcal{G}_\tau$ and $\hat{\mathbf{u}}$ denotes an extension of $\mathbf{u}$ from $\Omega^\delta$ into $\Omega$, see e.g. [42]. If $\Gamma_D = \emptyset$, then

$$V_{\delta} = \{\mathbf{u} \in H^1(\Omega^{\delta})^d : \int_{\Omega^{\delta}} \mathbf{u} \, dx = 0, \int_{\Omega} (\partial_{x_j} \hat{\mathbf{u}}_i - \partial_{x_i} \hat{\mathbf{u}}_j) \, dx = 0, \text{ for } i \neq j, i, j = 1, \dots, d\}.$$

Using (10) for $P^\delta$ and $\text{div}(J_g^\delta \mathbf{F}_{g,\delta}^{-T}) = 0$ in $\Omega^\delta$ we can rewrite equations in (7) as

$$\begin{align} -\operatorname{div}(J_g^\delta (\mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) + P_1(t, x)\mathbf{I}) \mathbf{F}_{g,\delta}^{-T}) &= -J_g^\delta \mathbf{F}_{g,\delta}^{-T} \nabla P_1(t, x) \quad \text{in } \Omega^\delta, t > 0, \tag{11} \\ J_g^\delta (\mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) + P_1(t, x)\mathbf{I}) \mathbf{F}_{g,\delta}^{-T} N &= -\delta J_g^\delta P_2^\delta(t, x) \mathbf{F}_{g,\delta}^{-T} N \quad \text{on } \Gamma^\delta, t > 0, \\ J_g^\delta \mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) \mathbf{F}_{g,\delta}^{-T} N &= J_g^\delta \mathbf{f}(t, x) |\mathbf{F}_{g,\delta}^{-T} N| \quad \text{on } \Gamma_N, t > 0, \\ \mathbf{u}^\delta \cdot \mathbf{N} = 0, \quad J_g^\delta \Pi_\tau (\mathbf{S}_\delta^T(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) \mathbf{F}_{g,\delta}^{-T} N) = 0 \quad \text{or} \quad \mathbf{u}^\delta = 0 \quad \text{on } \Gamma_D, t > 0, \end{align}$$

where $\Gamma_N = \partial\Omega \setminus \Gamma_D$ and $P_2^\delta(t, x) = P_2(t, x, x/\delta)$.

In the analysis and numerical implementation of (8) and (11) we shall consider weak solutions of the model equations. We shall use the notation $\langle\phi, \psi\rangle_A = \int_A \phi\psi \, dx$, for $\phi \in L^p(A)$, $\psi \in L^q(A)$, and $\langle\phi, \psi\rangle_{\partial A} = \int_{\partial A} \phi\psi \, d\gamma$, for $\phi \in L^p(\partial A)$, $\psi \in L^q(\partial A)$, with $1 < p, q < \infty$, $1/p + 1/q = 1$, and a bounded Lipschitz domain A.

**Definition 2.2.** A weak solution of (8) and (11) are $\mathbf{u}^\delta \in L^2(0, T; V_\delta)$ and $\mathbf{F}_{g,\delta} \in W^{1,\infty}(0, T; L^q(\Omega^\delta))^{d \times d}$, for any $q \in (1, \infty)$, with $\mathbf{F}_{g,\delta} \in L^\infty((0, T) \times \Omega^\delta)^{d \times d}$, satisfying

$$\begin{align} \left\langle J_g^\delta \mathbb{E}^\delta(x) \operatorname{sym}(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1}) + P_1(t, x) \mathbf{I}, \nabla \varphi \mathbf{F}_{g,\delta}^{-1} \right\rangle_{\Omega^\delta} + \delta \left\langle J_g^\delta P_2^\delta(t, x) \mathbf{F}_{g,\delta}^{-T} N, \varphi \right\rangle_{\Gamma^\delta} \tag{12} \\ - \left\langle J_g^\delta P_1(t, x) \mathbf{F}_{g,\delta}^{-T} N, \varphi \right\rangle_{\partial \Omega} = \left\langle J_g^\delta \mathbb{E}^\delta(x) [\mathbf{I} - \operatorname{sym}(\mathbf{F}_{g,\delta}^{-1})], \nabla \varphi \mathbf{F}_{g,\delta}^{-1} \right\rangle_{\Omega^\delta} \nonumber \\ - \left\langle J_g^\delta \mathbf{F}_{g,\delta}^{-T} \nabla P_1(t, x), \varphi \right\rangle_{\Omega^\delta} + \left\langle J_g^\delta \mathbf{f}(t, x) |\mathbf{F}_{g,\delta}^{-T} N|, \varphi \right\rangle_{\Gamma_N} \nonumber \end{align}$$

for $\varphi \in L^2(0, T; V_\delta)$ and a.a. $t \in (0, T)$, and $\mathbf{F}_{g,\delta}$ satisfies (8) a.e. in $(0, T) \times \Omega^\delta$.

First we shall prove a version of the Korn inequality, where the symmetric gradient includes the growth tensor.

**Lemma 2.3.** For $\mathbf{u} \in V_\delta$ and a tensor $\mathbf{F} \in L^\infty(\Omega)$ such that on each $\delta(Y + \xi) \cap \Omega$, for $\xi \in \Xi$, $\mathbf{F}$ is constant, $\det(\mathbf{F}) \ge 1$ and eigenvalues $\lambda_j(\mathbf{F}) \ge 1$ for $j = 1, \dots, d$, we have the following estimate

$$\|\mathbf{u}\|_{L^2(\Omega^\delta)}^2 + \|\nabla \mathbf{u}\|_{L^2(\Omega^\delta)}^2 \le C \|\operatorname{sym}(\nabla \mathbf{u} \mathbf{F}^{-1})\|_{L^2(\Omega^\delta)}^2, \quad (13)$$

with constant $C = C(\mathbf{F})$ independent of $\delta$.

**Proof.** Assumptions on the microscopic structure of the plant tissues ensure the existence of an extension $\hat{\mathbf{u}}$ of $\mathbf{u}$ from $\Omega^\delta$ to $\Omega$ with

$$\begin{align*} \|\hat{\mathbf{u}}\|_{H^1(\Omega)} &\le C \|\mathbf{u}\|_{H^1(\Omega^\delta)}, & \|\hat{\mathbf{u}}\|_{L^2(\Omega)} + \|\operatorname{sym}(\nabla \hat{\mathbf{u}})\|_{L^2(\Omega)} &\le C (\|\mathbf{u}\|_{L^2(\Omega^\delta)} + \|\operatorname{sym}(\nabla \mathbf{u})\|_{L^2(\Omega^\delta)}), \\ \|\nabla \hat{\mathbf{u}}\|_{L^2(\Omega)} &\le C \|\nabla \mathbf{u}\|_{L^2(\Omega^\delta)}, & \|\operatorname{sym}(\nabla \hat{\mathbf{u}})\|_{L^2(\Omega)} &\le C \|\operatorname{sym}(\nabla \mathbf{u})\|_{L^2(\Omega^\delta)}, \end{align*}$$
where constant $C$ is independent of $\delta$, see e.g. [42]. In the following we shall identify $\mathbf{u}$ with its extension $\hat{\mathbf{u}}$. Using that $\mathbf{F}$ and hence $\mathbf{F}^{-T}$ are constant on each $\delta(Y + \xi) \cap \Omega$, for $\xi \in \Xi$, and the properties of $\Omega^\delta$, together with an extension of $\mathbf{u}$ from $\delta\mathbf{F}(\xi)Y_w$ into $\delta\mathbf{F}(\xi)Y$, constructed as in e.g. [42, Lemma 4.1], we obtain

$$
\begin{align*}
\int_{\Omega^\delta} |\operatorname{sym}(\nabla \mathbf{u} \mathbf{F}^{-1})|^2 dx &= \sum_{\xi \in \Xi_\delta} \int_{\delta(Y_w + \xi)} |\operatorname{sym}(\nabla \mathbf{u} \mathbf{F}^{-1})|^2 dx + \sum_{l=1}^{L} \int_{\Lambda_\delta^l} |\operatorname{sym}(\nabla \mathbf{u} \mathbf{F}^{-1})|^2 dx \\
&= \sum_{\xi \in \Xi_\delta} \det(\mathbf{F}^{-1}(\xi)) \int_{\delta\mathbf{F}(\xi)(Y_w + \xi)} |\operatorname{sym}(\nabla \mathbf{u})|^2 dx + \sum_{l=1}^{L} \det(\mathbf{F}^{-1}(\xi_l)) \int_{\mathbf{F}(\xi_l)\Lambda_\delta^l} |\operatorname{sym}(\nabla \mathbf{u})|^2 dx \\
&\geq \sum_{\xi \in \Xi_\delta} \det(\mathbf{F}^{-1}(\xi)) C(\mathbf{F}(\xi)) \int_{\delta\mathbf{F}(\xi)(Y + \xi)} |\operatorname{sym}(\nabla \mathbf{u})|^2 dx + \sum_{l=1}^{L} \det(\mathbf{F}^{-1}(\xi_l)) \int_{\mathbf{F}(\xi_l)\Lambda_\delta^l} |\operatorname{sym}(\nabla \mathbf{u})|^2 dx,
\end{align*}
$$

where the dependence of the constant $C$ on $\mathbf{F}(\xi)$, for $\xi \in \Xi_\delta$, arises from the application of the Korn and Poincaré inequalities when constructing an extension from $\mathbf{F}(\xi)Y_w$ into $\mathbf{F}(\xi)Y$, see [42, Lemma 4.1] for more details, and $\xi_l \in \Lambda_\delta^l$ for $l = 1, \dots, L$. Using the uniform boundedness of $\mathbf{F}$, together with the fact that $\det(\mathbf{F}) \ge 1$ in $(0, T) \times \Omega$ and eigenvalues $\lambda_j(\mathbf{F}) \ge 1$, with $j = 1, \dots, d$, and applying the Korn inequality, see e.g. [10, 18, 42], yield

$$
\int_{\Omega^\delta} |\operatorname{sym}(\nabla \mathbf{u} \mathbf{F}^{-1})|^2 dx \geq C \int_{\Omega_F} |\operatorname{sym}(\nabla \mathbf{u})|^2 dx \geq C \int_{\Omega} |\operatorname{sym}(\nabla \mathbf{u})|^2 dx \\
\geq C \int_{\Omega} (|\mathbf{u}|^2 + |\nabla \mathbf{u}|^2) dx \geq C \int_{\Omega^\delta} (|\mathbf{u}|^2 + |\nabla \mathbf{u}|^2) dx,
$$

for $\mathbf{u} \in V_\delta$, where $\Omega_F = (\cup_{\xi \in \Xi_\delta} \delta\mathbf{F}(\xi)(Y + \xi)) \cup (\cup_{l=1}^L \mathbf{F}(\xi_l)\Lambda_\delta^l)$. This implies the result stated in the lemma. $\square$

**Remark 2.1.** Some results on the generalisation of the Korn inequality can be found in [38, 46]. Notice that in [38] $C^1$-regularity of $\mathbf{F}$ is required and in [46] the continuity of $\mathbf{F}$ or $\mathbf{u} = 0$ on $\partial\Omega$ are assumed. In the proof of Lemma 2.3 we use the fact that $\mathbf{F}$ is uniformly bounded and piece-wise constant and the eigenvalues of $\mathbf{F}$ satisfy $\lambda_j(\mathbf{F}) \ge 1$, for $j = 1, \dots, d$, ensuring that $\Omega$ is a subdomain of the transformed domain $\Omega_F$.

Using the inequality (13) and applying the Banach fixed point theorem we prove the well-posedness result for microscopic model (8) and (11).

**Theorem 2.4.** Under Assumptions 2.1 there exists a unique weak solution of (8), (11) satisfying

$$
(14) \qquad \|\mathbf{u}^\delta\|_{L^\infty(0,T;V_\delta)} + \|\mathbf{F}_{g,\delta}\|_{L^\infty((0,T)\times\Omega^\delta)} + \|\partial_t\mathbf{F}_{g,\delta}\|_{L^\infty((0,T)\times\Omega^\delta)} + \|\mathbf{F}_{g,\delta}\|_{W^{1,\infty}(0,T;L^q(\Omega^\delta))} \le C,
$$
with a constant $C$ independent of $\delta$ and any $q \in (1, \infty)$.

*Proof.* We shall apply a fixed point argument to show existence of a solution of (8), (11). Consider

$$
\mathcal{W} = \{ \mathbf{F} \in L^\infty((0, T) \times \Omega^\delta)^{d \times d} : \mathbf{F} \in W^{1,\infty}(0, T; L^q(\Omega^\delta))^{d \times d}, \|\mathbf{F}\|_{L^\infty((0,T) \times \Omega^\delta)} \le \hat{\kappa}, \det(\mathbf{F}) \ge 1, \\
\lambda_j(\mathbf{F}) \ge 1 \text{ in } [0, T] \times \Omega^\delta, j = 1, \dots, d, \text{ and piece-wise constant in each } \delta(Y_w + \xi) \cap \Omega \text{ for } \xi \in \Xi \}
$$

with $\hat{\kappa} \ge \exp(dMT)$. For a given $\tilde{\mathbf{F}}_{g,\delta} \in \mathcal{W}$, taking $\mathbf{u}^\delta$ as a test function in (12) and using assumptions on $\mathbb{E}$, $P_1$, $P_2$, and $\mathbf{f}$, together with the fact that $1 \le \det(\tilde{\mathbf{F}}_{g,\delta}) \le C$, with a constant $C$ independent of $\delta$, we obtain

$$
(15) \qquad \|\operatorname{sym}(\nabla \mathbf{u}^\delta \tilde{\mathbf{F}}_{g,\delta}^{-1})\|_{L^\infty(0,T;L^2(\Omega^\delta))}^2 \le C(1+\zeta\|\mathbf{u}^\delta\|_{L^\infty(0,T;V_\delta)}^2),
$$

for any fixed $\zeta > 0$. Here we also used the trace estimate

$$
\delta \|v\|^2_{L^2(\Gamma^\delta)} \le C (\|v\|^2_{L^2(\Omega^\delta)} + \delta^2 \|\nabla v\|^2_{L^2(\Omega^\delta)}) \quad \text{for } v \in H^1(\Omega^\delta),
$$

for some constant $C$ independent of $\delta$, see e.g. [27]. Notice that since $\tilde{\mathbf{F}}_{g,\delta}$ is constant on each $\delta(Y_w + \xi) \cap \Omega$ it can be trivially extended by the constant into each $\delta(Y + \xi) \cap \Omega$, for $\xi \in \Xi$. Using Lemma 2.3 we obtain

$$
(16) \qquad \|\mathbf{u}^\delta\|_{L^\infty(0,T;V_\delta)} \le C \|\operatorname{sym}(\nabla \mathbf{u}^\delta \tilde{\mathbf{F}}_{g,\delta}^{-1})\|_{L^\infty(0,T;L^2(\Omega^\delta))},
$$

and then choosing in (15) sufficiently small $\zeta > 0$ yields

$$
(17) \qquad \|\mathbf{u}^\delta\|_{L^\infty(0,T;V_\delta)} \le C,
$$
where the constant C does not depend on $\delta$. Assumptions on $\mathbb{E}$, $P_1, P_2$, $\mathbf{f}$ and $\tilde{\mathbf{F}}_{g,\delta}$, together with the estimate in Lemma 2.3, ensure that

$$ F(\varphi) = \langle J_g^\delta P_1(t, x) \tilde{\mathbf{F}}_{g,\delta}^{-T} N, \varphi \rangle_{\partial\Omega} - \langle P_1(t, x) \mathbf{I}, \nabla\varphi \tilde{\mathbf{F}}_{g,\delta}^{-1} \rangle_{\Omega^\delta} - \delta \langle J_g^\delta P_2(t, x) \tilde{\mathbf{F}}_{g,\delta}^{-T} N, \varphi \rangle_{\Gamma^\delta} \\ + \langle J_g^\delta \mathbb{E}^\delta(x) [\mathbf{I} - \text{sym}(\tilde{\mathbf{F}}_{g,\delta}^{-1})], \nabla\varphi \tilde{\mathbf{F}}_{g,\delta}^{-1} \rangle_{\Omega^\delta} - \langle J_g^\delta \tilde{\mathbf{F}}_{g,\delta}^{-T} \nabla P_1(t, x), \varphi \rangle_{\Omega^\delta} + \langle J_g^\delta \mathbf{f}(t, x) |\tilde{\mathbf{F}}_{g,\delta}^{-T} N|, \varphi \rangle_{\Gamma_N} $$

defines a bounded linear functional on $V_\delta$ and

$$ B(\mathbf{u}^\delta, \varphi) = \langle J_g^\delta \mathbb{E}^\delta(x) \text{sym}(\nabla \mathbf{u}^\delta \tilde{\mathbf{F}}_{g,\delta}^{-1}), \text{sym}(\nabla \varphi \tilde{\mathbf{F}}_{g,\delta}^{-1}) \rangle_{\Omega^\delta} $$

is a coercive bilinear form on $V_\delta \times V_\delta$, for $t \in (0, T)$ and a given $\tilde{\mathbf{F}}_{g,\delta} \in W$, where $J_g^\delta = \det(\tilde{\mathbf{F}}_{g,\delta})$. Thus the Lax-Milgram theorem, see e.g. [21], yields existence of a unique solution of problem (11) in $L^\infty(0, T; V_\delta)$ for a given $\tilde{\mathbf{F}}_{g,\delta} \in W$ and each fixed $\delta > 0$. The boundedness of function $G^\delta$ also ensures existence of a solution of (8) with $\tilde{\mathbf{F}}_{g,\delta}$ instead of $\mathbf{F}_{g,\delta}$ in $G^\delta$.

To show existence of a unique solution of the full model (8) and (11) we need to show a contraction property of the map $\mathcal{K}: W \to W$, where $\mathbf{F}_{g,\delta} = \mathcal{K}(\tilde{\mathbf{F}}_{g,\delta})$ is a solution of problem (11) and (8), with $\tilde{\mathbf{F}}_{g,\delta}$ instead of $\mathbf{F}_{g,\delta}$ in equation (11) and in function $G^\delta$ in (8).

From equations for $\mathbf{F}_{g,\delta}$, using properties of $G^\delta$, we obtain

$$ \det(\mathbf{F}_{g,\delta}(t, x)) = \det(\mathbf{F}_{g,\delta}(0, x)) \det\left(\exp\left(\int_0^t G^\delta(x, \nabla \mathbf{u}^\delta, \tilde{\mathbf{F}}_{g,\delta}^{-1}) ds\right)\right) = \exp\left(\text{tr}\left(\int_0^t G^\delta(x, \nabla \mathbf{u}^\delta, \tilde{\mathbf{F}}_{g,\delta}^{-1}) ds\right)\right) \ge 1, $$ 
for $(t, x) \in [0, T] \times \overline{\Omega^\delta}$, since diagonal elements of $G^\delta$ are nonnegative. Assumptions on $G^\delta$ imply that $\mathbf{F}_{g,\delta}$ is piece-wise constant in each $\delta(Y_w + \xi) \cap \Omega$ for $\xi \in \Xi$, has eigenvalues greater than or equal to 1, and 
$$ (18) \qquad \|\mathbf{F}_{g,\delta}\|_{L^\infty((0,T)\times\Omega^\delta)} + \|\partial_t \mathbf{F}_{g,\delta}\|_{L^\infty((0,T)\times\Omega^\delta)} + \|\mathbf{F}_{g,\delta}\|_{W^{1,\infty}(0,T;L^q(\Omega^\delta))} \le C, $$

with a constant C independent of $\delta$. Since $J_g^\delta = \det(\mathbf{F}_{g,\delta}) \ge 1$ for $(t, x) \in [0, T] \times \overline{\Omega^\delta}$, we also have

$$ \|\mathbf{F}_{g,\delta}^{-1}\|_{L^\infty((0,T)\times\Omega^\delta)} + \|\partial_t \mathbf{F}_{g,\delta}^{-1}\|_{L^\infty((0,T)\times\Omega^\delta)} + \|\mathbf{F}_{g,\delta}^{-1}\|_{W^{1,\infty}(0,T;L^q(\Omega^\delta))} \le C. $$

Multiplying the difference of (8) for $\mathbf{F}_{g,\delta}^1$ and $\mathbf{F}_{g,\delta}^2$ by $\mathbf{F}_{g,\delta}^1 - \mathbf{F}_{g,\delta}^2$, integrating over time variable, as well as using the boundedness of $\tilde{\mathbf{F}}_{g,\delta}$ and the Lipschitz continuity and boundedness of $G^\delta$, and applying the Gronwall inequality imply

$$ (19) \qquad \begin{aligned} \|\mathbf{F}_{g,\delta}^1 - \mathbf{F}_{g,\delta}^2\|_{L^\infty((0,T)\times\Omega^\delta)}^2 &\le TC \|G^\delta(x, \nabla \mathbf{u}_1^\delta, (\tilde{\mathbf{F}}_{g,\delta}^1)^{-1}) - G^\delta(x, \nabla \mathbf{u}_2^\delta, (\tilde{\mathbf{F}}_{g,\delta}^2)^{-1})\|_{L^\infty((0,T)\times\Omega^\delta)}^2 \\ &\le T[C_\delta \|\text{sym}(\nabla \mathbf{u}_1^\delta (\tilde{\mathbf{F}}_{g,\delta}^1)^{-1}) - \text{sym}(\nabla \mathbf{u}_2^\delta (\tilde{\mathbf{F}}_{g,\delta}^2)^{-1})\|_{L^\infty(0,T;L^2(\Omega^\delta))}^2 + C \|\tilde{\mathbf{F}}_{g,\delta}^1 - \tilde{\mathbf{F}}_{g,\delta}^2\|_{L^\infty((0,T)\times\Omega^\delta)}^2], \end{aligned} $$

for $\tilde{\mathbf{F}}_{g,\delta}^1, \tilde{\mathbf{F}}_{g,\delta}^2 \in W$, where constant C depends on $\|\tilde{\mathbf{F}}_{g,\delta}^j\|_{L^\infty((0,T)\times\Omega^\delta)}$, for $j = 1, 2$, which is bounded by $\hat{\kappa}$, and $C_\delta$ depends on $\|\tilde{\mathbf{F}}_{g,\delta}^j\|_{L^\infty((0,T)\times\Omega^\delta)}$ and $\delta$. In the derivation of (19) we used the following estimate

$$ \begin{aligned} & \|G^\delta(x, \nabla \mathbf{u}_1^\delta, (\tilde{\mathbf{F}}_{g,\delta}^1)^{-1}) - G^\delta(x, \nabla \mathbf{u}_2^\delta, (\tilde{\mathbf{F}}_{g,\delta}^2)^{-1})\|_{L^\infty(\Omega^\delta)}^2 \\ & \le C \left\| \int_{\delta([x/\delta]_Y + Y_w) \cap \Omega} \left| \boldsymbol{\varepsilon}^{\text{el}}(\nabla \mathbf{u}_1^\delta, (\tilde{\mathbf{F}}_{g,\delta}^1)^{-1}) - \boldsymbol{\varepsilon}^{\text{el}}(\nabla \mathbf{u}_2^\delta, (\tilde{\mathbf{F}}_{g,\delta}^2)^{-1}) \right| d\tilde{x} \right\|_{L^\infty(\Omega^\delta)}^2 \\ & \le C_1 \delta^{-d} \left\| \text{sym}(\nabla \mathbf{u}_1^\delta (\tilde{\mathbf{F}}_{g,\delta}^1)^{-1}) - \text{sym}(\nabla \mathbf{u}_2^\delta (\tilde{\mathbf{F}}_{g,\delta}^2)^{-1}) \right\|_{L^2(\Omega^\delta)}^2 + C_2 \|\tilde{\mathbf{F}}_{g,\delta}^1 - \tilde{\mathbf{F}}_{g,\delta}^2\|_{L^\infty(\Omega^\delta)}^2, \end{aligned} $$

for a.a. $t \in (0, T)$. From (11), using the uniform ellipticity of $\mathbb{E}$ and estimate (17), we obtain

$$ (20) \qquad \|\text{sym}(\nabla \mathbf{u}_1^\delta (\tilde{\mathbf{F}}_{g,\delta}^1)^{-1}) - \text{sym}(\nabla \mathbf{u}_2^\delta (\tilde{\mathbf{F}}_{g,\delta}^2)^{-1})\|_{L^\infty(0,T;L^2(\Omega^\delta))}^2 \le C \|\tilde{\mathbf{F}}_{g,\delta}^1 - \tilde{\mathbf{F}}_{g,\delta}^2\|_{L^\infty((0,T)\times\Omega^\delta)}^2. $$

In the derivation of (20) we used the following estimate

$$ \begin{aligned} & \|\nabla(\mathbf{u}_1^\delta - \mathbf{u}_2^\delta)[(\tilde{\mathbf{F}}_{g,\delta}^1)^{-1} - (\tilde{\mathbf{F}}_{g,\delta}^2)^{-1}]\|_{L^2(\Omega^\delta)} \le C_\zeta \|(\tilde{\mathbf{F}}_{g,\delta}^1)^{-1}\|_{L^\infty(\Omega^\delta)}^2 \|(\tilde{\mathbf{F}}_{g,\delta}^2)^{-1}\|_{L^\infty(\Omega^\delta)}^2 \|\tilde{\mathbf{F}}_{g,\delta}^1 - \tilde{\mathbf{F}}_{g,\delta}^2\|_{L^\infty(\Omega^\delta)}^2 \\ & + \zeta \|\nabla(\mathbf{u}_1^\delta - \mathbf{u}_2^\delta)\|_{L^2(\Omega^\delta)}^2 \le \zeta \|\text{sym}(\nabla(\mathbf{u}_1^\delta - \mathbf{u}_2^\delta)(\tilde{\mathbf{F}}_{g,\delta}^1)^{-1})\|_{L^2(\Omega^\delta)}^2 + C_1 \|\tilde{\mathbf{F}}_{g,\delta}^1 - \tilde{\mathbf{F}}_{g,\delta}^2\|_{L^\infty(\Omega^\delta)}^2 \\ & \le \zeta \|\text{sym}(\nabla\mathbf{u}_1^\delta(\tilde{\mathbf{F}}_{g,\delta}^1)^{-1}) - \text{sym}(\nabla\mathbf{u}_2^\delta(\tilde{\mathbf{F}}_{g,\delta}^2)^{-1})\|_{L^2(\Omega^\delta)}^2 + C(1 + \|\nabla\mathbf{u}_2^\delta\|_{L^2(\Omega^\delta)}^2) \|\tilde{\mathbf{F}}_{g,\delta}^1 - \tilde{\mathbf{F}}_{g,\delta}^2\|_{L^\infty(\Omega^\delta)}^2, \end{aligned} $$

for any fixed $\zeta > 0$. Here we applied (13) with $(\tilde{\mathbf{F}}_{g,\delta}^1)^{-1}$ and $\mathbf{u}_1^\delta(t) - \mathbf{u}_2^\delta(t) \in V_\delta$.

Combining estimates (19) and (20) and considering $T$ sufficiently small we obtain that $\mathcal{K}$ is a contraction. Then applying the Banach fixed point theorem yields existence of the unique solution for (8) and (11). Since
the choice of $T$ depends only on the parameters in the system and on $\delta$, and does not depend on the solution, we can iterate over time intervals to obtain the global existence and uniqueness result for the microscopic model (8) and (11) for any fixed $\delta > 0$. $\square$

Next we prove convergence results for sequences $\{\mathbf{u}^\delta\}$ and $\{\mathbf{F}_{g,\delta}\}$ as $\delta \to 0$.

**Lemma 2.5.** Under Assumptions 2.1, for sequence of solutions $\{\mathbf{u}^\delta\}$ and $\{\mathbf{F}_{g,\delta}\}$ of microscopic model (8) and (11), up to a subsequence, we have the following convergence results

$$
\begin{align}
\mathbf{u}^{\delta} &\to \mathbf{u} \qquad \text{weakly-* in } L^{\infty}(0, T; V), \tag{21} \\
\nabla \mathbf{u}^{\delta} &\to \nabla \mathbf{u} + \nabla_y \mathbf{u}_1 \qquad \text{two-scale, } \mathbf{u}_1 \in L^2((0, T) \times \Omega; H_{\text{per}}^1(Y)/\mathbb{R}), \\
\mathbf{F}_{g,\delta} &\to \mathbf{F}_g \qquad \text{weakly-* in } L^{\infty}((0, T) \times \Omega), \\
\mathbf{F}_{g,\delta} &\to \mathbf{F}_g \qquad \text{weakly-* in } W^{1,\infty}(0, T; L^q(\Omega)), \quad \text{for any } q \in (1, \infty),
\end{align}
$$

where we identify $\mathbf{u}^\delta$ with its extension from $\Omega^\delta$ into $\Omega$ and for $\mathbf{F}_{g,\delta}$ consider a trivial (constant) extension from $\delta(Y_w + \xi)$ into $\delta(Y + \xi)$ for $\xi \in \Xi_\delta$. The space $V$ is defined in the same way as $V_\delta$, with $\Omega^\delta$ replaced by $\Omega$. Assuming additionally that

$$ (22) \qquad \| \mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta}) \|_{L^{\infty}((0,T) \times \Omega; L^2(Y_w))} \le C, $$

with a constant $C$ independent of $\delta$, we have

$$
\begin{align}
& \mathbf{F}_{g,\delta} \to \mathbf{F}_g \qquad \text{strongly in } L^2((0,T) \times \Omega), \tag{23} \\
(23) \quad & \int_{Y_w} \mathcal{T}^{\delta}(\operatorname{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) dy \to \int_{Y_w} \operatorname{sym}((\nabla \mathbf{u} + \nabla_y \mathbf{u}_1) \mathbf{F}_g^{-1}) dy \qquad \text{strongly in } L^2((0,T) \times \Omega), \nonumber \\
& \int_{Y_w} \mathcal{T}^{\delta}(\mathbb{E}^{\delta}(x) \operatorname{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) dy \to \int_{Y_w} \mathbb{E}(\cdot, y) \operatorname{sym}((\nabla \mathbf{u} + \nabla_y \mathbf{u}_1) \mathbf{F}_g^{-1}) dy \qquad \text{strongly in } L^2((0,T) \times \Omega), \nonumber
\end{align}
$$

where $\mathcal{T}^\delta$ is the periodic unfolding operator defined below.

*Proof.* The first four convergences follow directly from estimates (14) and compactness results for the weak-* and two-scale convergences, see e.g. [3, 41] or Appendix for the definition and properties of the two-scale convergence.

To show the strong convergence of $\mathbf{F}_{g,\delta}$ we shall use the periodic unfolding operator $\mathcal{T}^\delta : L^p((0,T) \times \Omega^\delta) \to L^p((0,T) \times \Omega \times Y_w)$, with $1 \le p \le \infty$, defined as

$$ (24) \qquad \mathcal{T}^{\delta}(v)(t, x, y) = \begin{cases} v(t, \delta[x/\delta]_Y + \delta y) & \text{for a.e. } (x, y) \in (\Omega \setminus \Lambda_\delta) \times Y_w, t \in (0, T), \\ 0 & \text{for a.e. } (x, y) \in \Lambda_\delta \times Y_w, t \in (0, T), \end{cases} $$

where $[x/\delta]_Y$ denotes the unique integer combination, such that $x/\delta - [x/\delta]_Y$ belongs to $Y_w$ for $x \in \Omega^\delta \setminus \Lambda_\delta$, see e.g. [11, 13] and Appendix for more details.

Applying the periodic unfolding operator $\mathcal{T}^\delta$ to (8) yields

$$ (25) \qquad \partial_t \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) = G^{\delta}(x, \mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta}), \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}^{-1})) \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \quad \text{in } (0,T) \times \Omega \times Y_w, $$

where

$$ \tilde{G}^{\delta}(x, \mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta}), \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}^{-1})) = \eta_{\sigma} \left[ \int_{Y_w} \sigma(x, y, \mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta}), \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}^{-1})) dy - \tau_{\sigma} \right]_{+}, $$

or

$$ \tilde{G}^{\delta}(x, \mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta}), \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}^{-1})) = \eta_{\varepsilon} \left[ \int_{Y_w} \varepsilon^{\text{el}}(\mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta}), \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}^{-1})) dy - \tau_{\varepsilon} \right]_{+}, $$

and $G^\delta$ is defined in terms of $\tilde{G}^\delta$ as in (9). Notice that $G^\delta$ is independent of $y$ and hence $\mathcal{T}^\delta(\mathbf{F}_{g,\delta})$ depends on $t \in (0, T)$ and $x \in \Omega$ and is constant in $y \in Y_w$.

We shall show that $\{\mathcal{T}^\delta(\mathbf{F}_{g,\delta})\}$ converges strongly in $L^2((0, T) \times \Omega \times Y_w)$ by applying the Fréchet-Kolmogorov compactness theorem, [37], and its modification proposed by Simon, [52]. The uniform boundedness of $G^\delta$ implies

$$ (26) \qquad \| \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \|_{L^{\infty}((0,T) \times \Omega \times Y_w)} + \| \partial_t \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \|_{L^{\infty}((0,T) \times \Omega \times Y_w)} \le C, $$

where the constant $C$ is independent of $\delta$.
Denote $\Omega_h = \{x \in \Omega : \text{dist}(\partial\Omega, x) > h\}$ for $h > 0$. Considering (8) for $x + \hat{h}$ and $x$, with $\hat{h} = \delta\xi$, for $\xi \in \Xi$, and $|\hat{h}| \le h$, applying the unfolding operator $\mathcal{T}^\delta$, multiplying the resulting equations by $\mathcal{T}^\delta(\mathbf{F}_{g,\delta}(x + \hat{h})) - \mathcal{T}^\delta(\mathbf{F}_{g,\delta}(x))$, integrating over $(0,t)$, taking supremum over $Y_w$, and then integrating over $\Omega_{3h}$, we obtain

$$ (27) \qquad \begin{aligned} & \|\mathcal{T}^\delta(\mathbf{F}_{g,\delta}(\cdot + \hat{h})) - \mathcal{T}^\delta(\mathbf{F}_{g,\delta})\|_{L^2(\Omega_{3h}; L^\infty(Y_w))}^2 \le C \int_0^t \left( \|\mathcal{T}^\delta(\mathbf{F}_{g,\delta}(\cdot + \hat{h})) - \mathcal{T}^\delta(\mathbf{F}_{g,\delta})\|_{L^2(\Omega_{3h}; L^\infty(Y_w))}^2 \right. \\ & \qquad \left. + \left\| \int_{Y_w} \left[ \boldsymbol{\sigma}(\cdot + \hat{h}, y, \mathcal{T}^\delta(\nabla \mathbf{u}^\delta(\cdot + \hat{h})), \mathcal{T}^\delta(\mathbf{F}_{g,\delta}^{-1}(\cdot + \hat{h}))) - \boldsymbol{\sigma}(\cdot, y, \mathcal{T}^\delta(\nabla \mathbf{u}^\delta), \mathcal{T}^\delta(\mathbf{F}_{g,\delta}^{-1})) \right] dy \right\|_{L^2(\Omega_{3h})}^2 \right) ds, \end{aligned} $$

or with $\varepsilon^{\text{el}}$ instead of $\boldsymbol{\sigma}$ if we consider the second case for the growth rate. Notice that for the ease of notations in the formulas here and below we often do not write explicitly the dependence of $\mathbf{F}_{g,\delta}$, $\mathbf{u}^\delta$, $P_1$ and $P_2$ on $t \in [0, T]$.

To estimate the second term on the right-hand side of (27), we consider equations (11) for $x$ and $x + \hat{h}$ and take $(\mathbf{u}^\delta(t, x + \hat{h}) - \mathbf{u}^\delta(t, x))\rho_h^2(x)$ as a test function, where $\rho_h \in C_0^1(\Omega)$ with $\rho_h(x) = 1$ in $\Omega_{3h}$ and $\rho_h(x) = 0$ in $\Omega \setminus \Omega_{2h}$. Then applying the periodic unfolding operator $\mathcal{T}^\delta$ and using (22), together with the uniform boundedness of $\mathbf{F}_{g,\delta}$ and assumptions on $\mathbb{E}$, $P_1$, and $P_2$, yield

$$ (28) \qquad \begin{aligned} & \left\langle \mathcal{T}^\delta(J_g^\delta \mathbf{F}_{g,\delta}^{-1} \mathbb{E}^\delta(x + \hat{h})) \operatorname{sym}(\mathcal{T}^\delta(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1}(x + \hat{h}))) - \mathcal{T}^\delta(J_g^\delta \mathbf{F}_{g,\delta}^{-1} \mathbb{E}^\delta(x)) \operatorname{sym}(\mathcal{T}^\delta(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1}(x))) \right. \\ & \qquad \left. + \mathcal{T}^\delta(J_g^\delta \mathbf{F}_{g,\delta}^{-1} P_1(x + \hat{h})) - \mathcal{T}^\delta(J_g^\delta \mathbf{F}_{g,\delta}^{-1} P_1(x)), [\mathcal{T}^\delta(\nabla \mathbf{u}^\delta(x + \hat{h})) - \mathcal{T}^\delta(\nabla \mathbf{u}^\delta(x))] \mathcal{T}^\delta(\rho_h^2) \right\rangle_{\Omega_{2h} \times Y_w} \\ & \qquad + \left\langle \mathcal{T}^\delta(J_g^\delta P_2^\delta \mathbf{F}_{g,\delta}^{-T} N(x + \hat{h})) - \mathcal{T}^\delta(J_g^\delta P_2^\delta \mathbf{F}_{g,\delta}^{-T} N(x)), [\mathcal{T}^\delta(\mathbf{u}^\delta(x + \hat{h})) - \mathcal{T}^\delta(\mathbf{u}^\delta(x))] \mathcal{T}^\delta(\rho_h^2) \right\rangle_{\Omega_{2h} \times \Gamma} \\ & = - \left\langle \mathcal{T}^\delta(J_g^\delta \mathbf{F}_{g,\delta}^{-T} \nabla P_1(x + \hat{h})) - \mathcal{T}^\delta(J_g^\delta \mathbf{F}_{g,\delta}^{-T} \nabla P_1(x)), [\mathcal{T}^\delta(\mathbf{u}^\delta(x + \hat{h})) - \mathcal{T}^\delta(\mathbf{u}^\delta(x))] \mathcal{T}^\delta(\rho_h^2) \right\rangle_{\Omega_{2h} \times Y_w} \\ & \qquad + \left\langle \mathcal{T}^\delta(J_g^\delta \mathbf{F}_{g,\delta}^{-1} \mathbb{E}^\delta(x + \hat{h})) [\mathbf{I} - \operatorname{sym}(\mathcal{T}^\delta(\mathbf{F}_{g,\delta}^{-1}(x + \hat{h})))] \right. \\ & \qquad \left. - \mathcal{T}^\delta(J_g^\delta \mathbf{F}_{g,\delta}^{-1} \mathbb{E}^\delta(x)) [\mathbf{I} - \operatorname{sym}(\mathcal{T}^\delta(\mathbf{F}_{g,\delta}^{-1}(x)))], [\mathcal{T}^\delta(\nabla \mathbf{u}^\delta(x + \hat{h})) - \mathcal{T}^\delta(\nabla \mathbf{u}^\delta(x))] \mathcal{T}^\delta(\rho_h^2) \right\rangle_{\Omega_{2h} \times Y_w} \\ & \qquad + \kappa(h), \end{aligned} $$

where $\kappa(h) \to 0$ as $h \to 0$. Here we used that $|\hat{h}| \le h$, that $\rho_h(x) = 0$ and $\rho_h(x \pm \hat{h}) = 0$ for $x \in \Lambda_\delta$, and the following estimate

$$ (29) \qquad \begin{aligned} & \|\mathcal{T}^\delta(\mathbf{u}^\delta(\cdot + \hat{h}) - \mathbf{u}^\delta)\mathcal{T}^\delta(\rho_h \nabla \rho_h)\|_{L^2(\Omega_{2h} \times Y_w)} \\ & \le C \frac{|\hat{h}|}{h} \|\mathcal{T}^\delta(\nabla \mathbf{u}^\delta)\|_{L^2((\Omega_{2h} \setminus \Omega_{3h}) \times Y_w)} \le C |\Omega_{2h} \setminus \Omega_{3h}|^{1/2} \|\mathcal{T}^\delta(\nabla \mathbf{u}^\delta)\|_{L^\infty((0,T) \times \Omega; L^2(Y_w))} \le \kappa(h), \end{aligned} $$

for a.a. $t \in (0, T)$. Then using (22), together with the uniform boundedness of $\mathbf{F}_{g,\delta}$, the assumptions on $\mathbb{E}$, $P_1$, and $P_2$, and that integrals over $(\Omega_{2h} \setminus \Omega_{3h}) \times Y_w$ can be estimated by $|\Omega_{2h} \setminus \Omega_{3h}|$, from (28) we obtain

$$ (30) \qquad \begin{aligned} & \|\operatorname{sym}(\mathcal{T}^\delta(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1}(\cdot + \hat{h})) - \mathcal{T}^\delta(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1}))\|_{L^2(\Omega_{3h} \times Y_w)}^2 \le \kappa(h) \\ & \qquad + C \|\mathcal{T}^\delta(\mathbf{F}_{g,\delta}(\cdot + \hat{h})) - \mathcal{T}^\delta(\mathbf{F}_{g,\delta})\|_{L^2(\Omega_{3h}; L^\infty(Y_w))}^2, \end{aligned} $$

for a.a. $t \in (0, T)$, where $\kappa(h) \to 0$ as $h \to 0$ and $|\hat{h}| \le h$. Similar as in the proof of Theorem 2.4, in the derivation of (30) we used the following estimate

$$ \begin{aligned} & \|\mathcal{T}^\delta(\nabla \mathbf{u}^\delta(\cdot + \hat{h}) - \nabla \mathbf{u}^\delta) \mathcal{T}^\delta(\rho_h) [\mathcal{T}^\delta(\mathbf{F}_{g,\delta}^{-1}(\cdot + \hat{h})) - \mathcal{T}^\delta(\mathbf{F}_{g,\delta}^{-1})] \|_{L^1(\Omega_{2h}; L^2(Y_w))} \\ & \le \tilde{\zeta} \|\nabla(\mathbf{u}^\delta(\cdot + \hat{h}) - \mathbf{u}^\delta)\rho_h\|_{L^2(\Omega^\delta)} + C_\zeta \|\mathcal{T}^\delta(\mathbf{F}_{g,\delta}^{-1}(\cdot + \hat{h})) - \mathcal{T}^\delta(\mathbf{F}_{g,\delta}^{-1})\|_{L^2(\Omega_{2h}; L^\infty(Y_w))}^2 \\ & \le \zeta \|\operatorname{sym}(\nabla(\mathbf{u}^\delta(\cdot + \hat{h}) - \mathbf{u}^\delta)\mathbf{F}_{g,\delta}^{-1})\rho_h\|_{L^2(\Omega^\delta)}^2 + C_\zeta \|\mathcal{T}^\delta(\mathbf{F}_{g,\delta}(\cdot + \hat{h})) - \mathcal{T}^\delta(\mathbf{F}_{g,\delta})\|_{L^2(\Omega_{2h}; L^\infty(Y_w))}^2 \\ & \le \zeta \|\operatorname{sym}(\mathcal{T}^\delta(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1}(\cdot + \hat{h})) - \mathcal{T}^\delta(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1}))\|_{L^2(\Omega_{2h} \times Y_w)}^2 \\ & \qquad + C \|\mathcal{T}^\delta(\mathbf{F}_{g,\delta}(\cdot + \hat{h})) - \mathcal{T}^\delta(\mathbf{F}_{g,\delta})\|_{L^2(\Omega_{2h}; L^\infty(Y_w))}^2 + \kappa(h), \end{aligned} $$
for any fixed $\zeta > 0$, a.a. $t \in (0, T)$, and $\kappa(h) \to 0$ as $h \to 0$. In the derivation of the last estimate we used an estimate similar to (29), the properties of the unfolding operator and estimate (13) applied to $(\mathbf{u}^{\delta}(\cdot + \hat{h}) - \mathbf{u}^{\delta})\rho_h$ and $\mathbf{F}_{g,\delta}^{-1}$. Notice that $(\mathbf{u}^{\delta}(\cdot + \hat{h}) - \mathbf{u}^{\delta})\rho_h = 0$ on $\partial\Omega$.

Using (30) in (27), together with the regularity of $\mathbb{E}$ with respect to the first variable, and applying the Grönwal inequality yields

$$ (31) \qquad \| \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}(\cdot, \cdot + \hat{h})) - \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \|_{L^{\infty}(0,T; L^2(\Omega_{3h}; L^{\infty}(Y_w)))}^2 \le \kappa(h), $$

and, using (30), also

$$ (32) \qquad \| \operatorname{sym}(\mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1}(\cdot, \cdot + \hat{h})) - \mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) \|_{L^{\infty}(0,T; L^2(\Omega_{3h} \times Y_w))}^2 \le \kappa(h), $$

where $\kappa(h) \to 0$ as $h \to 0$ and $|\hat{h}| \le h$. In a similar way, using the uniform boundedness of $\partial_t \mathbf{F}_{g,\delta}$, and hence also of $\partial_t \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta})$, and uniform continuity, with respect to the time variable, of $P_1$, $\nabla P_1$, $P_2$, and $\mathbf{f}$, we obtain

$$ (33) \qquad \| \operatorname{sym}(\mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1}(\cdot, \cdot + \tilde{\tau}, \cdot) - \mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) \|_{L^2((0,T-\tilde{\tau}) \times \Omega \times Y_w)}^2 \le \tilde{\kappa}(\tilde{\tau}), \quad \text{with } \tilde{\kappa}(\tilde{\tau}) \to 0 \text{ as } \tilde{\tau} \to 0. $$

From the definition of the periodic unfolding operator, for any $\tilde{h} \in \mathbb{R}^d$, with $|\tilde{h}| < h$, we have

$$ \begin{aligned} \| \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}(\cdot, \cdot + \tilde{h}, \cdot) - \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \|_{L^2(\Omega_{3h,T} \times Y_w)}^2 &\le |Y| \sum_{k \in \{0,1\}^d} \| \mathbf{F}_{g,\delta}(\cdot, \cdot + \delta(k_b + [\tilde{h}/\delta]_Y)) - \mathbf{F}_{g,\delta} \|_{L^2(\Omega_{3h,T}^\delta)}^2 \\ &\le \sum_{k \in \{0,1\}^d} \| \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}(\cdot, \cdot + \hat{h})) - \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \|_{L^2(\Omega_{3h,T} \times Y_w)}^2 \le \sum_{k \in \{0,1\}^d} \| \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}(\cdot, \cdot + \hat{h})) - \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \|_{L^2(\Omega_{3h,T}; L^{\infty}(Y_w))}^2, \end{aligned} $$

where $\hat{h} = \delta(k_b + [\tilde{h}/\delta]_Y)$, with $k_b = \sum_{j=1}^d k_j b_j$ and $|\hat{h}| \le h$ for sufficiently small $\delta$, and $\Omega_{3h,T} = (0,T) \times \Omega_{3h}$. Then for any $h > 0$ there exists $\delta_0$, such that $|\hat{h}| = |\delta(k_b + [\tilde{h}/\delta]_Y)| \le h$ for all $\delta \le \delta_0$. Considering (31) for such $\hat{h}$ yields

$$ (34) \qquad \| \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}(\cdot, \cdot + \tilde{h}, \cdot) - \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \|_{L^2((0,T) \times \Omega_{3h} \times Y_w)}^2 \le \kappa(h), $$

where $\kappa(h) \to 0$ as $h \to 0$, and hence as $|\tilde{h}| \to 0$, for all $\delta \le \delta_0$. For $\delta > \delta_0$ we have a finite number of members of the sequence and for each of them the continuity of the $L^2$-norm of $L^2$-functions, see e.g. [6], ensures (34) for an appropriate $\tilde{h}$, with $|\tilde{h}| < h$. Then from the finite number of such $\tilde{h}$ considering the smallest one implies the property (34) for all $\delta > 0$. Similarly we obtain

$$ (35) \qquad \| \operatorname{sym}(\mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1}(\cdot, \cdot + \tilde{h}, \cdot)) - \mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) \|_{L^2((0,T) \times \Omega_{3h} \times Y_w)}^2 \le \kappa(h), $$

for $\tilde{h} \in \mathbb{R}^d$, with $|\tilde{h}| < h$, and all $\delta > 0$, where $\kappa(h) \to 0$ as $h \to 0$. Thus using the uniform boundedness of $\mathbf{F}_{g,\delta}$, assumption (22), and the fact that $\mathcal{T}^{\delta}(\mathbf{F}_{g,\delta})$ is constant in $y \in Y_w$, and applying the Fréchet-Kolmogorov and Simon compactness theorems, see [52], we obtain the strong convergence of $\mathcal{T}^{\delta}(\mathbf{F}_{g,\delta})$ to $\tilde{\mathbf{F}}_g$ in $L^2((0,T) \times \Omega \times Y_w)$ and of $\int_{Y_w} \operatorname{sym}(\mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) dy$ in $L^2((0,T) \times \Omega)$. The properties of $\mathcal{T}^{\delta}$ and that $\tilde{\mathbf{F}}_g$ and $\mathcal{T}^{\delta}(\mathbf{F}_{g,\delta})$ are constant in $y \in Y_w$ imply

$$ \begin{aligned} \int_{\Omega_T} \mathbf{F}_g \psi \, dxdt &= \lim_{\delta \to 0} \int_{\Omega_T} \mathbf{F}_{g,\delta} \psi \, dxdt = \frac{1}{|Y|} \lim_{\delta \to 0} \int_{\Omega_T} \int_Y \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \mathcal{T}^{\delta}(\psi) \, dydxdt \\ &= \frac{1}{|Y_w|} \lim_{\delta \to 0} \int_{\Omega_T} \int_{Y_w} \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) \mathcal{T}^{\delta}(\psi) \, dydxdt = \frac{1}{|Y_w|} \int_{\Omega_T} \int_{Y_w} \tilde{\mathbf{F}}_g \psi \, dydxdt = \int_{\Omega_T} \tilde{\mathbf{F}}_g \psi \, dxdt, \end{aligned} $$

for all $\psi \in L^2((0,T) \times \Omega)$. Hence $\tilde{\mathbf{F}}_g = \mathbf{F}_g$ in $(0,T) \times \Omega$. Similarly we obtain

$$ \begin{aligned} \| \mathbf{F}_{g,\delta} - \mathbf{F}_g \|_{L^2((0,T) \times \Omega)}^2 &\le |Y|^{-1} \| \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) - \mathcal{T}^{\delta}(\mathbf{F}_g) \|_{L^2((0,T) \times \Omega \times Y)}^2 + C|\Lambda_{\delta}| \\ &\le |Y_w|^{-1} \| \mathcal{T}^{\delta}(\mathbf{F}_{g,\delta}) - \mathbf{F}_g \|_{L^2((0,T) \times \Omega \times Y_w)}^2 + |Y|^{-1} \| \mathcal{T}^{\delta}(\mathbf{F}_g) - \mathbf{F}_g \|_{L^2((0,T) \times \Omega \times Y)}^2 + C|\Lambda_{\delta}| \to 0, \quad \text{as } \delta \to 0, \end{aligned} $$

which ensures the strong convergence of $\mathbf{F}_{g,\delta}$ to $\mathbf{F}_g$ in $L^2((0,T) \times \Omega)$. Strong convergence of $\mathbf{F}_{g,\delta}$, together with the uniform boundedness of $\mathbf{F}_{g,\delta}$ and $\det(\mathbf{F}_{g,\delta}) \ge 1$, implies also strong convergence of $\mathbf{F}_{g,\delta}^{-1}$ in $L^2((0,T) \times \Omega)$. Then using the equivalence between two-scale convergence of a sequence and weak convergence of the corresponding unfolded sequence we obtain $\int_{Y_w} \operatorname{sym}(\mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) dy \to \int_{Y_w} \operatorname{sym}((\nabla \mathbf{u} + \nabla_y \mathbf{u}_1) \mathbf{F}_g^{-1}) dy$ strongly in $L^2((0,T) \times \Omega)$.

Since the elasticity tensor $\mathbb{E}$ is independent of the time variable, estimate (33) implies

$$ (36) \qquad \| \mathcal{T}^{\delta}(\mathbb{E}^{\delta}\operatorname{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) (\cdot + \tilde{\tau}, \cdot, \cdot) - \mathcal{T}^{\delta}(\mathbb{E}^{\delta}\operatorname{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) \|_{L^2((0,T-\tilde{\tau}) \times \Omega \times Y_w)}^2 \le \tilde{\kappa}(\tilde{\tau}), $$
with $\tilde{\kappa}(\tilde{\tau}) \to 0$ as $\tilde{\tau} \to 0$. From (35), together with the regularity of $\mathbb{E}$ with respect to the first variable and $Y$-periodicity with respect to the second variable, we obtain

$$ \begin{aligned} & \| \mathcal{T}^\delta(\mathbb{E}^\delta)(\cdot + \tilde{h}, \cdot) \operatorname{sym}(\mathcal{T}^\delta(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1}))(\cdot, \cdot + \tilde{h}, \cdot) - \mathcal{T}^\delta(\mathbb{E}^\delta) \operatorname{sym}(\mathcal{T}^\delta(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1})) \|_{L^2((0,T) \times \Omega_{3h} \times Y_w)}^2 \\ & \leq \|\mathbb{E}\|_{L^\infty(\Omega \times Y_w)}^2 \kappa(h) + \sum_{k \in \{0,1\}^d} \|\mathbb{E}(\cdot + \hat{h}, \cdot) - \mathbb{E}\|_{L^\infty(\Omega_{3h} \times Y_w)}^2 \|\operatorname{sym}(\nabla \mathbf{u}^\delta \mathbf{F}_{g,\delta}^{-1})\|_{L^2((0,T) \times \Omega^\delta)}^2 \leq \kappa_1(h), \end{aligned} $$

where $\hat{h} = \delta(k_b + [\tilde{h}/\delta]_Y)$ for $\tilde{h} \in \mathbb{R}^d$ with $|\tilde{h}| < h$ and $\kappa_1(h) \to 0$ as $h \to 0$. Then using the two-scale convergence of $\nabla \mathbf{u}^\delta$ to $\nabla \mathbf{u} + \nabla_y \mathbf{u}_1$, the strong convergence of $\mathbf{F}_{g,\delta}^{-1}$, and the strong two-scale convergence of $\mathbb{E}^\delta$ to $\mathbb{E}$, together with the Fréchet-Kolmogorov compactness theorem, yields the last convergence in (23). $\square$

### 3. DERIVATION OF MACROSCOPIC EQUATIONS FOR MICROSCOPIC MODEL (8) AND (11)

We shall use both the formal asymptotic expansion and two-scale convergence methods to derive macroscopic equations for (8),(11).

#### 3.1. Formal asymptotic expansion.
First we present the formal derivation of the macroscopic equations using the asymptotic expansion of $\mathbf{u}^\delta$ and $\mathbf{F}_{g,\delta}^{-1}$ in powers of $\delta$. The convergence results (21) ensure that the limit functions (zero-order terms) $\mathbf{u} = \mathbf{u}(t, x)$ and $\mathbf{F}_g = \mathbf{F}_g(t, x)$ are independent of the microscopic variable $y \in Y_w$. Hence the formal asymptotic expansion ansatz reads

$$ \begin{aligned} (37) \qquad & \mathbf{u}^\delta(t, x) = \mathbf{u}(t, x) + \delta\mathbf{u}_1(t, x, x/\delta) + \delta^2\mathbf{u}_2(t, x, x/\delta) + \dots, \\ & \mathbf{F}_{g,\delta}^{-1}(t, x) = \mathbf{F}_g^{-1}(t, x) + \delta\mathbf{F}_{g,1}^{-1}(t, x, x/\delta) + \delta^2\mathbf{F}_{g,2}^{-1}(t, x, x/\delta) + \dots, \end{aligned} $$

for $(t, x) \in (0, T) \times \Omega^\delta$ and $\mathbf{u}_j(t, x, \cdot)$ and $\mathbf{F}_{g,j}^{-1}(t, x, \cdot)$ are $Y$-periodic, for $j = 1, 2, \dots$. Substituting (37) into microscopic equations (11) yields

$$ (38) \qquad -(\operatorname{div}_x + \frac{1}{\delta}\operatorname{div}_y) \left( J_g^\delta \left[ \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}} \left( \nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1 + \delta(\nabla_x \mathbf{u}_1 + \nabla_y \mathbf{u}_2) \dots, \mathbf{F}_g^{-1} + \delta \mathbf{F}_{g,1}^{-1} \dots \right) \right. \right. \\ \qquad \left. \left. + P_1(t, x) \mathbf{I} \right] (\mathbf{F}_g^{-T} + \delta \mathbf{F}_{g,1}^{-T} \dots) \right) = -J_g^\delta (\mathbf{F}_g^{-T} + \delta \mathbf{F}_{g,1}^{-T} \dots) \nabla_x P_1(t, x) $$

in $(0, T) \times \Omega \times Y_w$, where $J_g^\delta = J_g^\delta((\mathbf{F}_g^{-1} + \delta \mathbf{F}_{g,1}^{-1})^{-1})$. For the boundary conditions on $(0, T) \times \Omega \times \Gamma$ we have

$$ (39) \qquad J_g^\delta \left[ \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}} \left( \nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1 + \delta(\nabla_x \mathbf{u}_1 + \nabla_y \mathbf{u}_2) \dots, \mathbf{F}_g^{-1} + \delta \mathbf{F}_{g,1}^{-1} \dots \right) \right. \\ \qquad \left. + P_1(t, x) \mathbf{I} \right] (\mathbf{F}_g^{-T} + \delta \mathbf{F}_{g,1}^{-T} \dots) N = -\delta J_g^\delta P_2(t, x, y) (\mathbf{F}_g^{-T} + \delta \mathbf{F}_{g,1}^{-T} \dots) N $$

and on $(0, T) \times \Gamma_N \times Y_w$ we obtain

$$ (40) \qquad J_g^\delta \left[ \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}} (\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1 + \delta(\nabla_x \mathbf{u}_1 + \nabla_y \mathbf{u}_2) \dots, \mathbf{F}_g^{-1} + \delta \mathbf{F}_{g,1}^{-1} \dots) + P_1(t, x) \mathbf{I} \right] (\mathbf{F}_g^{-T} + \delta \mathbf{F}_{g,1}^{-T} \dots) N \\ = J_g^\delta \left[ P_1(t, x) (\mathbf{F}_g^{-T} + \delta \mathbf{F}_{g,1}^{-T} \dots) N + \mathbf{f}(t, x) |(\mathbf{F}_g^{-T} + \delta \mathbf{F}_{g,1}^{-T} \dots) N| \right]. $$

For the growth equation it holds

$$ (41) \quad \partial_t (\mathbf{F}_g^{-1} + \delta \mathbf{F}_{g,1}^{-1} \dots)^{-1} = G^\delta (x, \nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1 + \delta(\nabla_x \mathbf{u}_1 + \nabla_y \mathbf{u}_2) \dots, \mathbf{F}_g^{-1} + \delta \mathbf{F}_{g,1}^{-1} \dots) (\mathbf{F}_g^{-1} + \delta \mathbf{F}_{g,1}^{-1} \dots)^{-1}. $$

Now we shall consider terms for different powers of $\delta$. For $O(\delta^{-1})$ in equation (38) and $O(1)$ in boundary condition (39) we obtain

$$ (42) \qquad \begin{aligned} -\operatorname{div}_y & \left[ J_g \left( \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}} (\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) + P_1(t, x) \mathbf{I} \right) \mathbf{F}_g^{-T} \right] = 0 && \text{in } (0, T) \times \Omega \times Y_w, \\ J_g & \left( \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}} (\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) + P_1(t, x) \mathbf{I} \right) \mathbf{F}_g^{-T} N = 0 && \text{on } (0, T) \times \Omega \times \Gamma, \\ \mathbf{u}_1 & && \text{Y-periodic,} \end{aligned} $$

where $J_g = J_g(\mathbf{F}_g)$. This is an elliptic problem for $\mathbf{u}_1$ in $Y_w$ for given $\mathbf{u}$ and $\mathbf{F}_g$.
Terms of O(1) in equations (38) are

$$ 
\begin{aligned}
& - \operatorname{div}_x \left[ J_g \left( \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) + P_1(t, x) \mathbf{I} \right) \mathbf{F}_g^{-T} \right] \\
& - \operatorname{div}_y \left[ J'_g \cdot \mathbf{F}_{g,1}^{-1} \left( \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) + P_1(t, x) \mathbf{I} \right) \mathbf{F}_g^{-T} \right]
\end{aligned}
$$

$$ (43) \quad - \operatorname{div}_y \left[ J_g \left( \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) + P_1(t, x) \mathbf{I} \right) \mathbf{F}_{g,1}^{-T} \right] \\ - \operatorname{div}_y \left[ J_g \mathbb{E}(x, y) \partial_2 \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) \mathbf{F}_{g,1}^{-1} \mathbf{F}_g^{-T} \right] \\ - \operatorname{div}_y \left[ J_g \mathbb{E}(x, y) \partial_1 \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) (\nabla_x \mathbf{u}_1 + \nabla_y \mathbf{u}_2) \mathbf{F}_g^{-T} \right] = -J_g \mathbf{F}_g^{-T} \nabla_x P_1(t, x) \quad \text{in } \Omega \times Y_w $$

and O(δ)-terms in boundary condition (39) are

$$ 
\begin{aligned}
& J'_g \cdot \mathbf{F}_{g,1}^{-1} \left( \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) + P_1(t, x) \mathbf{I} \right) \mathbf{F}_g^{-T} N \\
& + J_g \left( \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) + P_1(t, x) \mathbf{I} \right) \mathbf{F}_{g,1}^{-T} N \\
& + J_g \mathbb{E}(x, y) \partial_1 \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) (\nabla_x \mathbf{u}_1 + \nabla_y \mathbf{u}_2) \mathbf{F}_g^{-T} N \\
& + J_g \mathbb{E}(x, y) \partial_2 \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) \mathbf{F}_{g,1}^{-1} \mathbf{F}_g^{-T} N = -J_g P_2(t, x, y) \mathbf{F}_g^{-T} N \quad \text{on } \Omega \times \Gamma,
\end{aligned}
$$

where $J'_g \cdot \mathbf{F}_{g,1}^{-1}$ denotes the sum of the components of the Hadamard product of two matrices $J'_g$ and $\mathbf{F}_{g,1}^{-1}$.

Equations (43), with boundary condition (44) and Y-periodicity of $\mathbf{u}_2$, define an elliptic problem for $\mathbf{u}_2$ in $Y_w$ for given $\mathbf{u}, \mathbf{u}_1, \mathbf{F}_g^{-1}$, and $\mathbf{F}_{g,1}^{-1}$. Applying the Fredholm alternative to ensure existence of a solution $\mathbf{u}_2$ of (43) and (44), see e.g. [21], and using Y-periodicity of $\mathbb{E}, \mathbf{F}_{g,1}^{-1}$, and $\mathbf{u}_j$, for $j = 1, 2$, imply macroscopic equations for $\mathbf{u}$:

$$ 
\begin{aligned}
& -\operatorname{div}_x \int_{Y_w} J_g \left( \mathbb{E}(x, y) \boldsymbol{\varepsilon}^{\mathrm{el}}(\nabla_x \mathbf{u} + \nabla_y \mathbf{u}_1, \mathbf{F}_g^{-1}) + P_1(t, x) \mathbf{I} \right) \mathbf{F}_g^{-T} dy \\
& = -|Y_w| J_g \mathbf{F}_g^{-T} \nabla_x P_1(t, x) - \int_{\Gamma} J_g P_2(t, x, y) \mathbf{F}_g^{-T} N d\gamma_y \quad \text{in } (0, T) \times \Omega.
\end{aligned}
$$

As next we need to determine $\mathbf{u}_1$. Using the expression for $\boldsymbol{\varepsilon}^{\mathrm{el}}$ and transformation of problem (42) to be defined on $Y_{w,g} = \mathbf{F}_g Y_w$ and $\Gamma_g = \mathbf{F}_g \Gamma$ for $(t, x) \in (0, T) \times \Omega$, we obtain

$$ 
\begin{aligned}
& \operatorname{div}_y \left( \mathbb{E}(x, \mathbf{F}_g^{-1}y) [\operatorname{sym}(\nabla_y \mathbf{u}_1) + \operatorname{sym}(\nabla_x \mathbf{u}\mathbf{F}_g^{-1}) + \operatorname{sym}(\mathbf{F}_g^{-1}) - \mathbf{I}] + P_1(t, x) \mathbf{I} \right) = 0 \quad \text{in } Y_{w,g}, \\
& (46) \quad \left( \mathbb{E}(x, \mathbf{F}_g^{-1}y) [\operatorname{sym}(\nabla_y \mathbf{u}_1) + \operatorname{sym}(\nabla_x \mathbf{u}\mathbf{F}_g^{-1}) + \operatorname{sym}(\mathbf{F}_g^{-1}) - \mathbf{I}] + P_1(t, x) \mathbf{I} \right) \nu = 0 \quad \text{on } \Gamma_g, \\
& \mathbf{u}_1 \quad Y_g - \text{periodic.}
\end{aligned}
$$

This is a linear elliptic problem for $\mathbf{u}_1$ which due to assumptions on $\mathbb{E}$ has, up to constant in $y$, a unique solution. Thus we can consider $\mathbf{u}_1$ in the form

$$ (47) \quad \mathbf{u}_1(t, x, y) = \sum_{i,j=1}^{d} \left[ (\nabla_x \mathbf{u} \mathbf{F}_g^{-1})_{ij}(t, x) + (\mathbf{F}_g^{-1}(t, x) - \mathbf{I})_{ij} \right] \mathbf{w}^{ij}(t, x, y) + P_1(t, x) \mathbf{v}(t, x, y) + \bar{\mathbf{u}}_1(t, x), $$

where $\mathbf{w}^{ij}$ and $\mathbf{v}$ are solutions of the 'unit cell' problems

$$ 
\begin{aligned}
& \operatorname{div}_y \left( \mathbb{E}(x, \mathbf{F}_g^{-1}y) [\operatorname{sym}(\nabla_y \mathbf{w}^{ij}) + \mathbf{b}_{ij}] \right) = 0 \quad \text{in } Y_{w,g}, \\
& (48) \quad \mathbb{E}(x, \mathbf{F}_g^{-1}y) [\operatorname{sym}(\nabla_y \mathbf{w}^{ij}) + \mathbf{b}_{ij}] \nu = 0 \quad \text{on } \Gamma_g, \quad \mathbf{w}^{ij} \quad Y_g - \text{periodic},
\end{aligned}
$$

for $i, j = 1, \dots, d$, where $\mathbf{b}_{ij} = \frac{1}{2}(\mathbf{e}_i \otimes \mathbf{e}_j + \mathbf{e}_j \otimes \mathbf{e}_i)$ and $\{\mathbf{e}_j\}_{j=1}^d$ is the standard basis in $\mathbb{R}^d$, and

$$ 
\begin{aligned}
& \operatorname{div}_y \left( \mathbb{E}(x, \mathbf{F}_g^{-1}y) \operatorname{sym}(\nabla_y \mathbf{v}) + \mathbf{I} \right) = 0 \quad \text{in } Y_{w,g}, \\
& (49) \quad \left( \mathbb{E}(x, \mathbf{F}_g^{-1}y) \operatorname{sym}(\nabla_y \mathbf{v}) + \mathbf{I} \right) \nu = 0 \quad \text{on } \Gamma_g, \quad \mathbf{v} \quad Y_g - \text{periodic.}
\end{aligned}
$$

Transforming back to the fixed domain $Y_w$ gives

$$ 
\begin{aligned}
& \operatorname{div}_y \left( J_g \mathbb{E}(x, y) (\operatorname{sym}(\nabla_y \mathbf{w}^{ij} \mathbf{F}_g^{-1}) + \mathbf{b}_{ij}) \mathbf{F}_g^{-T} \right) = 0 \quad \text{in } Y_w, \\
& J_g \mathbb{E}(x, y) (\operatorname{sym}(\nabla_y \mathbf{w}^{ij} \mathbf{F}_g^{-1}) + \mathbf{b}_{ij}) \mathbf{F}_g^{-T} N = 0 \quad \text{on } \Gamma, \quad \mathbf{w}^{ij} \quad Y - \text{periodic,}
\end{aligned}
$$
$$
\text{and}
$$
$$
(51) \qquad \operatorname{div}_y \left( J_g (\mathbb{E}(x, y)\operatorname{sym}(\nabla_y \mathbf{v}\mathbf{F}_g^{-1}) + \mathbf{I}) \mathbf{F}_g^{-T} \right) = 0 \quad \text{in } Y_w,
$$
$$
J_g (\mathbb{E}(x, y)\operatorname{sym}(\nabla_y \mathbf{v}\mathbf{F}_g^{-1}) + \mathbf{I}) \mathbf{F}_g^{-T} N = 0 \quad \text{on } \Gamma, \quad \mathbf{v}, \mathbf{Y} \text{ - periodic.}
$$

Using the solutions of the ‘unit cells’ problems, the macroscopic (homogenized) elasticity tensor E<sub>hom</sub> is defined
as

$$
(52) \quad \begin{aligned}
\mathbb{E}_{\text{hom},ijkl}(t,x) &= \frac{1}{|Y_g|} \int_{Y_{w,g}} \left( \mathbb{E}_{ijkl}(x, \mathbf{F}_g^{-1}y) + \left[ \mathbb{E}(x, \mathbf{F}_g^{-1}y) \operatorname{sym}(\nabla_y \mathbf{w}^{ij}) \right]_{kl} \right) dy \\
&= \frac{1}{|Y|} \int_{Y_w} \left( \mathbb{E}_{ijkl}(x,y) + \left[ \mathbb{E}(x,y) \operatorname{sym}(\nabla_y \mathbf{w}^{ij} \mathbf{F}_g^{-1}) \right]_{kl} \right) dy \quad \text{for } (t,x) \in (0,T) \times \Omega,
\end{aligned}
$$

and matrix $K_{\text{hom}}$ is given by

$$
(53) \quad K_{\text{hom}}(t,x) = \frac{1}{|Y_g|} \int_{Y_{w,g}} \mathbb{E}(x, \mathbf{F}_g^{-1}y) \operatorname{sym}(\nabla_y \mathbf{v}) dy = \frac{1}{|Y|} \int_{Y_w} \mathbb{E}(x, y) \operatorname{sym}(\nabla_y \mathbf{v} \mathbf{F}_g^{-1}) dy \quad \text{in } (0,T) \times \Omega.
$$

Using the structure of **u**<sub>1</sub> in (45) and in O(1)-terms in (40), together with the expressions for E<sub>hom</sub> and K<sub>hom</sub>, yields the macroscopic problem (54) formulated below. Considering O(1)-terms in (41) and using again the structure of **u**<sub>1</sub> imply macroscopic equations (55) for the growth tensor **F**<sub>g</sub>.

**3.2. Rigorous derivation of macroscopic model.** Now, using the two-scale convergence method, see e.g. [3, 14, 41], and the properties of the periodic unfolding operator, see e.g. [13], we rigorously derive the macroscopic model for the microscopic problem (8) and (11).

**Theorem 3.1.** Under assumptions in Lemma 2.5, up to a subsequence, solutions {$\mathbf{u}^\delta$} and {$\mathbf{F}_{g,\delta}$} of the mi-
croscopic problem (8) and (11) converge, as $\delta \to 0$, to solution $\mathbf{u} \in L^\infty(0, T; V)$ and $\mathbf{F}_g \in W^{1,\infty}(0, T; L^q(\Omega)) \cap L^\infty((0,T) \times \Omega)$, for any $q \in (1, \infty)$, of the macroscopic problem

$$
\begin{equation}
\begin{split}
\operatorname{div}_x \left( J_g \boldsymbol{\sigma}_{\text{hom}} \mathbf{F}_g^{-T} \right) &= J_g \mathbf{F}_g^{-T} \frac{1}{|Y|} \int_{\Gamma} P_2(t, x, y) N d\gamma_y \quad \text{in } (0, T) \times \Omega, \\
(54) \qquad J_g \boldsymbol{\sigma}_{\text{hom}} \mathbf{F}_g^{-T} N &= \left(1 - \frac{|Y_w|}{|Y|}\right) J_g P_1(t, x) \mathbf{F}_g^{-T} N + J_g |\mathbf{F}_g^{-T} N| \mathbf{f}(t, x) \quad \text{on } (0, T) \times \Gamma_N, \\
\mathbf{u} \cdot N = 0, \quad J_g \Pi_\tau (\boldsymbol{\sigma}_{\text{hom}} \mathbf{F}_g^{-T} N) &= \left(1 - \frac{|Y_w|}{|Y|}\right) J_g P_1(t, x) \Pi_\tau (\mathbf{F}_g^{-T} N), \quad \text{or } \mathbf{u} = 0 \quad \text{on } (0, T) \times \Gamma_D,
\end{split}
\end{equation}
$$

$$
\begin{align*}
& \text{and} \\
(55) \qquad & \partial_t \mathbf{F}_g = G(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \mathbf{F}_g && \text{in } (0, T) \times \Omega, \\
& \mathbf{F}_g(0) = \mathbf{I} && \text{in } \Omega,
\end{align*}
$$

$$
\begin{align*}
& \textit{with} \\
(56) \qquad & \tilde{G}(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) = \eta_{\sigma} [(|Y|/|Y_w|) \boldsymbol{\sigma}_{\text{hom}}(t, x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) - \boldsymbol{\tau}_{\sigma}]_+ && \text{or} \\
(57) \qquad & \tilde{G}(x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) = \eta_{\varepsilon} [(|Y|/|Y_w|) \boldsymbol{\varepsilon}_{\text{hom}}(t, x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) - \boldsymbol{\tau}_{\varepsilon}]_+, \nonumber
\end{align*}
$$

$$
\begin{align*}
& \textit{where} \\
(58) \qquad \boldsymbol{\sigma}_{\text{hom}}(t, x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) &= \mathbb{E}_{\text{hom}}(t, x) \boldsymbol{\varepsilon}^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g^{-1}) + K_{\text{hom}}(t, x) P_1(t, x) \\
& \qquad = \frac{|Y_w|}{|Y|} \boldsymbol{\varepsilon}^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g^{-1}) + \sum_{i,j=1}^{d} \boldsymbol{\varepsilon}^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g^{-1})_{ij} \frac{1}{|Y|} \int_{Y_w} \operatorname{sym}(\nabla_y \mathbf{w}^{ij} \mathbf{F}_g^{-1}) dy \\
(59) \qquad & \qquad + P_1(t, x) \frac{1}{|Y|} \int_{Y_w} \operatorname{sym}(\nabla_y \mathbf{v} \mathbf{F}_g^{-1}) dy,
\end{align*}
$$

with $\boldsymbol{\varepsilon}^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g^{-1}) = \text{sym}(\nabla \mathbf{u} \mathbf{F}_g^{-1}) + \text{sym}(\mathbf{F}_g^{-1}) - \mathbf{I}$, and $G$ is defined by $\tilde{G}$ as in (6). The macroscopic tensors $\mathbb{E}_{\text{hom}}$ and $K_{\text{hom}}$ are defined by (52) and (53), where $\mathbf{w}^{ij}$ and $\mathbf{v}$ are solutions of the 'unit cell' problems (50) and (51). The space $V$ is defined in the same way as $V_{\delta}$ with $\Omega^{\delta}$ replaced by $\Omega$.

Assuming, additionally,

$$
(60) \qquad \| \nabla \mathbf{u} \|_{L^{\infty}((0,T) \times \Omega)} \leq C,
$$

solution of problem (54), (55) is unique and the whole sequence of solutions of microscopic model (8) and (11)
converges to the solution of the macroscopic problem.
Proof. Considering $\varphi(t, x) = \psi(t, x) + \delta\psi_1(t, x, x/\delta)$, with $\psi \in C^1([0, T]; C^1(\Omega) \cap V)$ and $\psi_1 \in C_0^1((0, T) \times \Omega; C_{\text{per}}^1(Y))$, as a test function in (12), taking the two-scale limit as $\delta \to 0$, and using the two-scale convergence of $\nabla u^\delta$ and the strong convergence of $\mathbf{F}_{g,\delta}^{-1}$, ensured by the boundedness of $\mathbf{F}_{g,\delta}^{-1}$ and the strong convergence of $\mathbf{F}_{g,\delta}$, see Lemma 2.5, we obtain

$$ (61) \qquad \int_{\Omega_T} \frac{1}{|Y|} \int_{Y_w} J_g \left[ \mathbb{E}(x, y) (\operatorname{sym}(\nabla \mathbf{u} \mathbf{F}_g^{-1} + \nabla_y \mathbf{u}_1 \mathbf{F}_g^{-1}) + \operatorname{sym}(\mathbf{F}_g^{-1}) - \mathbf{I}) \right. \\ \qquad \qquad \left. + P_1(t, x) \mathbf{I} \right] (\nabla \psi + \nabla_y \psi_1) \mathbf{F}_g^{-1} dydxdt = - \int_{\Omega_T} \frac{1}{|Y|} \int_{Y_w} J_g \mathbf{F}_g^{-T} \nabla P_1(t, x) \psi dydxdt \\ \qquad \qquad + \int_{(\partial\Omega)_T} J_g P_1(t, x) \mathbf{F}_g^{-T} N \psi d\gamma dt + \int_{\Gamma_{N,T}} J_g \mathbf{f}(t, x) |\mathbf{F}_g^{-T} N| \psi d\gamma dt \\ \qquad \qquad \left. - \int_{\Omega_T} \frac{1}{|Y|} \int_{\Gamma} J_g P_2(t, x, y) \mathbf{F}_g^{-T} N \psi d\gamma_y dxdt, \right] $$

where $\Omega_T = (0, T) \times \Omega$, $(\partial\Omega)_T = (0, T) \times \partial\Omega$, and $\Gamma_{N,T} = (0, T) \times \Gamma_N$. Considering now $\psi = 0$ implies the problem for $\mathbf{u}_1$:

$$ (62) \qquad \int_{\Omega_T} \frac{1}{|Y|} \int_{Y_w} J_g \left( \mathbb{E}(x, y) \left[ \operatorname{sym}(\nabla_y \mathbf{u}_1 \mathbf{F}_g^{-1}) + \operatorname{sym}(\nabla \mathbf{u} \mathbf{F}_g^{-1}) + \operatorname{sym}(\mathbf{F}_g^{-1}) - \mathbf{I} \right] \right. \\ \qquad \qquad \left. + P_1(t, x) \mathbf{I} \right) \nabla_y \psi_1 \mathbf{F}_g^{-1} dydxdt = 0. $$

From (62), considering for $\mathbf{u}_1$ the ansatz (47) we obtain the ‘unit cell’ problems (48) and (49) and the formulas for the macroscopic coefficient $\mathbb{E}_{\text{hom}}$ and $K_{\text{hom}}$ in (52) and (53), respectively.

Consider $\psi_1 = 0$ in (61) yields

$$ (63) \qquad \int_{\Omega_T} J_g \left[ \frac{1}{|Y|} \int_{Y_w} \mathbb{E}(x, y) (\operatorname{sym}(\nabla \mathbf{u} \mathbf{F}_g^{-1} + \nabla_y \mathbf{u}_1 \mathbf{F}_g^{-1}) + \operatorname{sym}(\mathbf{F}_g^{-1}) - \mathbf{I}) dy + \frac{|Y_w|}{|Y|} P_1(t, x) \mathbf{I} \right] \nabla \psi \mathbf{F}_g^{-1} dxdt \\ = - \int_{\Omega_T} \frac{|Y_w|}{|Y|} J_g \mathbf{F}_g^{-T} \nabla P_1(t, x) \psi dxdt + \int_{(\partial\Omega)_T} J_g P_1(t, x) \mathbf{F}_g^{-T} N \psi d\gamma dt \\ + \int_{\Gamma_{N,T}} J_g \mathbf{f}(t, x) |\mathbf{F}_g^{-T} N| \psi d\gamma dt - \int_{\Omega_T} J_g \mathbf{F}_g^{-T} \frac{1}{|Y|} \int_{\Gamma} P_2(t, x, y) N d\gamma_y \psi dxdt. $$

Using the structure of $\mathbf{u}_1$ and expressions for $\mathbb{E}_{\text{hom}}$ and $K_{\text{hom}}$, we obtain

$$ (64) \qquad \int_{\Omega_T} J_g \left[ \mathbb{E}_{\text{hom}}(t, x) \left[ \operatorname{sym}(\nabla \mathbf{u} \mathbf{F}_g^{-1}) + \operatorname{sym}(\mathbf{F}_g^{-1}) - \mathbf{I} \right] + \left( K_{\text{hom}}(t, x) + \frac{|Y_w|}{|Y|} \mathbf{I} - \mathbf{I} \right) P_1(t, x) \right] \nabla \psi \mathbf{F}_g^{-1} dxdt \\ = \int_{\Omega_T} J_g \left( 1 - \frac{|Y_w|}{|Y|} \right) \mathbf{F}_g^{-T} \nabla P_1(t, x) \psi dxdt + \int_{\Gamma_{N,T}} J_g |\mathbf{F}_g^{-T} N| \mathbf{f}(t, x) \psi d\gamma dt \\ - \int_{\Omega_T} J_g \mathbf{F}_g^{-T} \frac{1}{|Y|} \int_{\Gamma} P_2(t, x, y) N d\gamma_y \psi dxdt, $$

which is the macroscopic problem (54) for $\mathbf{u}$ in the weak form. Here we used $\text{div}(J_g \mathbf{F}_g^{-T}) = 0$ in $(0, T) \times \Omega$ to rewrite the boundary term involving $P_1$ as an integral over $\Omega$ of $\nabla P_1$. Using the same calculations, the first integral on the right hand side in (64), combined with the last two terms on the left hand side, can be rewritten as the integral over $(0, T) \times \partial\Omega$.

To pass to the limit in the equation (8), we first determine the weak limit of $\int_{\delta([x/\delta]_Y+Y_w)} \sigma(\tilde{x}, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) d\tilde{x}$ and $\int_{\delta([x/\delta]_Y+Y_w)} \varepsilon^{\text{el}}(\nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) d\tilde{x}$. Using the two-scale convergence of $\nabla \mathbf{u}^\delta$ and strong convergence of $\mathbf{F}_{g,\delta}^{-1}$, together with the relation between the two-scale convergence of a sequence and weak convergence of the unfolded sequence, yields

$$ \begin{aligned} \lim_{\delta \to 0} \int_0^T \int_{\Omega} \frac{1}{|Y_w|} \int_{\delta([x/\delta]_Y+Y_w)} \sigma(\tilde{x}, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1}) d\tilde{x} \, \phi \, dxdt &= \lim_{\delta \to 0} \int_0^T \int_{\Omega} \frac{1}{|Y_w|} \int_{Y_w} \mathcal{T}^\delta(\sigma(x, \nabla \mathbf{u}^\delta, \mathbf{F}_{g,\delta}^{-1})) dy \, \phi \, dxdt \\ &= \int_0^T \int_{\Omega} \frac{1}{|Y_w|} \int_{Y_w} \mathbb{E}(x, y) [\operatorname{sym}(\nabla \mathbf{u} \mathbf{F}_g^{-1}) + \operatorname{sym}(\nabla_y \mathbf{u}_1 \mathbf{F}_g^{-1}) + \operatorname{sym}(\mathbf{F}_g^{-1}) - \mathbf{I}] dy \, \phi \, dxdt \\ &= \int_0^T \int_{\Omega} \frac{|Y|}{|Y_w|} \sigma_{\text{hom}}(t, x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \phi \, dxdt, \end{aligned} $$
for $\phi \in C_0((0,T) \times \Omega)$ and $\sigma_{\text{hom}}(t, x, \nabla \mathbf{u}, \mathbf{F}_g^{-1})$ given by (58). Similarly we obtain the weak convergence for the strain

$$ \lim_{\delta \to 0} \int_0^T \int_{\Omega} \frac{1}{|Y_w|} \int_{\delta([x/\delta]_Y + Y_w)} \varepsilon^{\text{el}}(\nabla \mathbf{u}^{\delta}, \mathbf{F}_{g,\delta}^{-1}) d\tilde{x} \, \phi \, dxdt = \int_0^T \int_{\Omega} \frac{|Y|}{|Y_w|} \varepsilon_{\text{hom}}(t, x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) \, \phi \, dxdt, $$

where the macroscopic strain $\varepsilon_{\text{hom}}(t, x, \nabla \mathbf{u}, \mathbf{F}_g^{-1})$ is defined by (59). Then the continuity of $G$ and the strong convergence of $\mathbf{F}_{g,\delta}$, $\int_{Y_w} \mathcal{T}^{\delta}(\text{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) dy$, and $\int_{Y_w} \mathcal{T}^{\delta}(\mathbb{E}^{\delta}(x)\text{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) dy$ yield macroscopic equation (55).

Considering a piece-wise constant approximation of $\mathbf{F}_g$ by $\mathbf{F}_{g,j} = \text{const}$ on $\Omega_j$, with $\Omega = \text{Int}(\cup_{j=1}^n \overline{\Omega}_j)$, $\Omega_i \cap \Omega_j = \emptyset$, for $i \neq j$, $J_{g,j} = \det(\mathbf{F}_{g,j}) \ge 1$ and eigenvalues $\lambda_k(\mathbf{F}_{g,j}) \ge 1$, for $k = 1, \dots, d$ and $j = 1, \dots, n$, and using the Korn inequality for $\mathbf{u}(t) \in V$, we can write

$$ \begin{aligned} \int_{\Omega} |\text{sym}(\nabla \mathbf{u} \mathbf{F}_g^{-1})|^2 dx &= \lim_{n \to \infty} \sum_{j=1}^n \int_{\Omega_j} |\text{sym}(\nabla \mathbf{u} \mathbf{F}_{g,j}^{-1})|^2 dx = \lim_{n \to \infty} \sum_{j=1}^n \int_{\mathbf{F}_{g,j} \Omega_j} J_{g,j}^{-1} |\text{sym}(\nabla \mathbf{u})|^2 dx \\ &\ge C \lim_{n \to \infty} \sum_{j=1}^n \int_{\Omega_j} |\text{sym}(\nabla \mathbf{u})|^2 dx = C \int_{\Omega} |\text{sym}(\nabla \mathbf{u})|^2 dx \ge C \|\mathbf{u}\|_V^2. \end{aligned} \tag{65} $$

Then using (60), in a similar way as in the proof of Lemma 2.5, for two solutions $\mathbf{u}_1, \mathbf{F}_{g,1}$ and $\mathbf{u}_2, \mathbf{F}_{g,2}$ of (54),(55) we obtain

$$ (66) \qquad \|\text{sym}(\nabla \mathbf{u}_1 \mathbf{F}_{g,1}^{-1}) - \text{sym}(\nabla \mathbf{u}_2 \mathbf{F}_{g,2}^{-1})\|_{L^2(\Omega)}^2 \le C \|\mathbf{F}_{g,1} - \mathbf{F}_{g,2}\|_{L^2(\Omega)}^2, $$

for a.a. $t \in (0, T)$. Considering then the equation for $\mathbf{F}_{g,1} - \mathbf{F}_{g,2}$ and using (66) yields $\|\mathbf{F}_{g,1}(t) - \mathbf{F}_{g,2}(t)\|_{L^2(\Omega)} = 0$ for a.a. $t \in (0, T)$, and hence $\mathbf{F}_{g,1} = \mathbf{F}_{g,2}$ a.e. in $(0, T) \times \Omega$. Using this in (66) implies $\text{sym}(\nabla(\mathbf{u}_1 - \mathbf{u}_2)\mathbf{F}_g^{-1}) = 0$ a.e. in $(0, T) \times \Omega$ and hence, due to (65), also $\mathbf{u}_1 = \mathbf{u}_2$ a.e. in $(0, T) \times \Omega$ and the uniqueness of solution of (54),(55). Thus we have that the whole sequence of solutions of the microscopic problem converges to the solution of the macroscopic equations. $\square$

**Remark 3.1.** (i) Using assumptions on the domain $\Omega$ and regularity of $\mathbb{E}$, $P_1$, $P_2$ and $\mathbf{f}$, estimates (22) and (60), assumed to be true in the rigorous derivation of the macroscopic model, may possibly be shown considering approaches similar to [1, 2, 9, 30, 51].

(ii) In the derivation of the strong convergence of $\mathbf{F}_{g,\delta}$ and $\int_{Y_w} \mathcal{T}^{\delta}(\text{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})) dy$ we assumed the uniform in $\delta$ boundedness (22) of the $L^2(Y_w)$-norm of $\mathcal{T}^{\delta}(\nabla \mathbf{u}^{\delta})$ in $(0, T) \times \Omega$.

It is also possible to derive macroscopic equations (61) for $\mathbf{u}$, by assuming first the strong convergence of $\mathbf{F}_{g,\delta}$ and then by showing the strong two-scale convergence of $\text{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})$ deduce the strong convergence of $\mathbf{F}_{g,\delta}$. To show the strong two-scale convergence of $\text{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})$, we consider $\nabla \mathbf{u}^{\delta}$ as a test function in (12) and take the limit as $\delta \to 0$. Then using the lower semi-continuity of the norm, together with the positivity of $J_g$ and properties of $\mathbb{E}$, yields

$$ \begin{aligned} & \frac{1}{|Y|} \left\langle J_g \mathbb{E}(x, y) \text{sym}((\nabla \mathbf{u} + \nabla_y \mathbf{u}_1) \mathbf{F}_g^{-1}), \text{sym}((\nabla \mathbf{u} + \nabla_y \mathbf{u}_1) \mathbf{F}_g^{-1}) \right\rangle_{\Omega_T, Y_w} \\ & \le \lim_{\delta \to 0} \left\langle J_g^{\delta} \mathbb{E}^{\delta}(x) \text{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1}), \text{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1}) \right\rangle_{\Omega_T^{\delta}} \\ & = -\frac{1}{|Y|} \left\langle J_g P_1(t, x) \mathbf{I}, (\nabla \mathbf{u} + \nabla_y \mathbf{u}_1) \mathbf{F}_g^{-1} \right\rangle_{\Omega_T, Y_w} - \frac{1}{|Y|} \left\langle J_g P_2(t, x, y) \mathbf{F}_g^{-T} N, \mathbf{u} \right\rangle_{\Omega_T, \Gamma} \\ & \quad + \left\langle J_g P_1(t, x) \mathbf{F}_g^{-T} N, \mathbf{u} \right\rangle_{(\partial \Omega)_T} + \frac{1}{|Y|} \left\langle J_g \mathbb{E}(x, y) [\mathbf{I} - \text{sym}(\mathbf{F}_g^{-1})], (\nabla \mathbf{u} + \nabla_y \mathbf{u}_1) \mathbf{F}_g^{-1} \right\rangle_{\Omega_T, Y_w} \\ & \quad - \frac{1}{|Y|} \left\langle J_g \mathbf{F}_g^{-T} \nabla P_1(t, x), \mathbf{u} \right\rangle_{\Omega_T, Y_w} + \left\langle J_g \mathbf{f}(x) |\mathbf{F}_g^{-T} N|, \mathbf{u} \right\rangle_{\Gamma_{N,T}} \\ & = \frac{1}{|Y|} \left\langle J_g \mathbb{E}(x, y) \text{sym}((\nabla \mathbf{u} + \nabla_y \mathbf{u}_1) \mathbf{F}_g^{-1}), \text{sym}((\nabla \mathbf{u} + \nabla_y \mathbf{u}_1) \mathbf{F}_g^{-1}) \right\rangle_{\Omega_T, Y_w}, \end{aligned} $$

where the last equality follows from (61) by considering $\psi = \mathbf{u}$ and $\psi_1 = \mathbf{u}_1$. This, together with the two-scale convergence of $\mathbb{E}^{\delta}(x)\text{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})$ and the strong convergence of $J_g^{\delta}$, ensured by the strong convergence of $\mathbf{F}_{g,\delta}$, implies the corresponding strong two-scale convergence of $\text{sym}(\nabla \mathbf{u}^{\delta} \mathbf{F}_{g,\delta}^{-1})$. Using now this result in equation (8) yields the strong convergence of $\mathbf{F}_{g,\delta}$.
4. NUMERICAL SIMULATIONS OF MICROSCOPIC AND MACROSCOPIC TWO-SCALE PROBLEMS

To numerically solve the microscopic problem (11) and (8) we consider a rectangular tissue of hexagonal
cells with geometrical parameters – wall length $l_1$, $l_2$, wall thickness $w$ and angle $\theta$ – as shown in the unit cell
scheme of Figure 1a-b. The tissue is composed of $N_x$ cells along the first axis and $N_x/2$ cell layers along the
second axis (see Figure 2d). The bottom cell layer is cut horizontally at mid-height, and we assume no normal
displacement and no tangential stress at the cut walls. As long as there is no symmetry breaking, this setting
amounts to simulating a square tissue made of $N_x(N_x - 1)$ cells. The parameters related to the microscopic
mechanical properties of the cell wall are its Young modulus $E$ and its Poisson ratio $\nu$. Regarding growth,
we denote for the numerical simulations simply by $\eta$ the extensibility $\eta_\sigma$ or $\eta_\epsilon$, and simply by $\tau$ the yield
threshold tensor $\tau_\sigma$ or $\tau_\epsilon$ in case of stress based and strain based growth respectively. Furthermore, in all our
simulations we consider an isotropic threshold tensor $\tau = \tau I$ with parameter $\tau$.

For the space-dependent turgor pressure $P_1(x)$ a piece-wise constant approximation $P^\delta$ of $P_1$ is defined as in (10), with an appropriate function $P_2$, where $P_2(x, \cdot)$ is $Y$- periodic. In the case of a linear function $P_1(x) = \alpha x_1 + \beta$ we have $P^\delta(x) = \beta$ for $x \in [0, \delta)$ and $P^\delta(x) = \beta + \alpha \xi_1 \delta$ for $x \in \delta(Y + \xi)$ with $\xi \in \Xi^\delta$ and $P_2(x, y) = -\alpha \hat{P}(y)$, where $\hat{P}(y) = y_1$ on $Y$ and $Y$-periodically extended to $\mathbb{R}^d$. Notice that for constant $P_1$ we have $P_2 = 0$.

In order to solve the microscopic problem at each time-step, we consider its variational formulation (12) and implement it using the Finite Element Method and the open-source finite element software FreeFEM [25]. The vector components of the displacement fields **u** and test function fields **v** are represented by $P1$ finite elements, and we use the sparsesolver solver of FreeFEM. The mesh is built so that a fixed number of triangle edges is imposed per unit length of domain boundary (in most of the simulations we use 25 triangle edges per cell wall length $l_2$).

The numerical simulations of the macroscopic two-scale problem are based on an algorithm that we schema-
tized in Figure 2b and detailed in Algorithm 1. Namely, at each time-step, we solve the macroscopic problem
(Algorithm 3) in order to compute tissue level growth and displacement. To compute the required tissue level
material properties, we locally solve the corresponding unit cell problems (Algorithm 2). In doing so, we use
two meshes (Figure 2a) – a fine mesh **m** for the numerical resolution of the macroscopic problem, and a coarse
mesh ℳ for the calculation of the homogenized (macroscopic) properties $E_{\text{hom}}$ and $K_{\text{hom}}$. For instance for the
simulations presented in Figures 2 and 3, we consider a tissue of width 7 and use 10 edges per unit length for
the fine mesh **m**, corresponding to 70 triangle edges along the width of the tissue. In contrast, the coarse mesh
**(71)** for the same tissue has edge length of 1, corresponding to 7 triangle edges along the width (see Figure 2a).

To determine $E_{\text{hom}}$ and $K_{\text{hom}}$ for given $\mathbf{F}_g(t, \hat{x})$, where $\hat{x}$ are the nodal points of the coarse mesh, we compute numerically solutions of the 'unit cell' problems (50) and (51), written in the weak form as

$$
(67) \quad \int_{Y_w} J_g \mathbb{E}(x, y) \left[ \text{sym}(\nabla_y \mathbf{w}^{ij} \mathbf{F}_g^{-1}) + \mathbf{b}_{ij} \right] \nabla_y \phi \mathbf{F}_g^{-1} dy = 0 \quad \text{for } \phi \in H_{\text{per}}^1(Y), \quad i, j = 1, \dots, d,
$$

$$
(68) \quad \int_{Y_w} J_g (\mathbb{E}(x, y) \operatorname{sym}(\nabla_y \mathbf{v} \mathbf{F}_g^{-1}) + \mathbf{I}) \nabla_y \phi \mathbf{F}_g^{-1} dy = 0 \quad \text{for } \phi \in H_{\text{per}}^1(Y).
$$

We call the solutions **w**<sup>ij</sup> and **v** of the ‘unit cell’ problems the elementary deformations (see Figure 1b),
as they form the building blocks for the effective material properties of the tissue. If we denote an arbitrary
displacement field by **w**, a generalised gradient and its product with the elastic tensor by

$$
d[\mathbf{w}](t, x, y) = \operatorname{sym}(\nabla_y \mathbf{w}(t, x, y) \mathbf{F}_g^{-1}(t, x)) \quad \text{and} \quad \mathbb{E}d[\mathbf{w}](t, x, y) = \mathbb{E}(x, y)\operatorname{sym}(\nabla_y \mathbf{w}(t, x, y) \mathbf{F}_g^{-1}(t, x)),
$$

and the mean over the ‘unit cell’ centered on the node with coordinates $x$ by

$$
\overline{d[\mathbf{w}]}(t, x) = \frac{1}{|Y|} \int_{Y_w} d[\mathbf{w}](t, x, y)dy \quad \text{and} \quad \overline{\mathbb{E}d[\mathbf{w}]}(t, x) = \frac{1}{|Y|} \int_{Y_w} \mathbb{E}d[\mathbf{w}](t, x, y)dy,
$$

then the effective material properties (52) and (53) are computed as

$$
(69) \quad \mathbb{E}_{\text{hom},ijkl} = \overline{E_{ijkl}} + (\overline{\mathbb{E}d[\mathbf{w}^{ij}]})_{kl}, \quad K_{\text{hom}} = \overline{\mathbb{E}d[\mathbf{v}]}. \qquad \text{by}
$$

On the same nodal points we also compute the components of the vector field

$$
(70) \qquad \mathbf{P}_2(t, x) = \frac{1}{|Y|} \int_{\Gamma} P_2(t, x, y) N d\gamma_y.
$$

Then we interpolate the values for $E_{\text{hom}}$, $K_{\text{hom}}$ and $\mathbf{P}_2$ to obtain the corresponding tensor and vector fields defined in the whole domain $\Omega$, so that these can be used in the numerical simulations of the macroscopic
equation (64). With **u** the solution of the macroscopic problem at time *t*, we compute the macroscopic strain $\varepsilon_{\text{hom}}$ and stress $\sigma_{\text{hom}}$ tensor fields,

$$ (71) \qquad \begin{aligned} \varepsilon_{\text{hom}}(t, x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) &= \bar{\mathbf{I}}(t, x)\varepsilon^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g^{-1}) + \sum_{i,j=1}^{d} \overline{d[\mathbf{w}^{ij}]}(t, x)\varepsilon^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g^{-1})_{ij} + \overline{d[\mathbf{v}]}(t, x)P_1(t, x), \\ \sigma_{\text{hom}}(t, x, \nabla \mathbf{u}, \mathbf{F}_g^{-1}) &= \mathbb{E}_{\text{hom}}(t, x)\varepsilon^{\text{el}}(\nabla \mathbf{u}, \mathbf{F}_g) + K_{\text{hom}}(t, x)P_1(t, x), \end{aligned} $$

and the corresponding growth rate at this time step according to the strain or stress hypothesis (57) or (56). Finally we use the Euler method to compute the growth tensor $\mathbf{F}_g$ for the next time-step according to equation (55).

**Algorithm 1: Coupled simulation**

**Inputs**: reference tissue geometry,
$P_1$, $E$, $\nu$, $\eta$, $\tau$ functions defined on the tissue domain
time step $dt$ and maximal time $t_{max}$

**Output**: $\mathbf{F}_g(t)$ field for all time points $t$ in $[0, t_{max}]$,
macroscopic displacement field $\mathbf{u}(t)$ for all time points $t$ in $[0, t_{max}]$.

1. Construct fine mesh **m** and coarse mesh $\mathfrak{M}$ of the tissue;

2. Initialize $t = 0$;

3. Initialize $\mathbf{F}_g = \mathbf{I}$ on $\mathfrak{M}$;

4. Evaluate the functions $E$, $\nu$ on each element of the coarse mesh $\mathfrak{M}$;

5. Define FE fields $P_1$, $\eta$ and $\tau$ on the fine mesh **m**;

6. **while** $t \le t_{max}$ **do**

7. &nbsp;&nbsp;&nbsp;&nbsp;**for** all elements of the coarse mesh $\mathfrak{M}$ **do**

8. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Solve** the ‘unit cell’ problems (Algorithm 2);

9. &nbsp;&nbsp;&nbsp;&nbsp;Reconstruct FE fields of the effective properties $\mathbb{E}_{\text{hom}}$, $K_{\text{hom}}$ and $\mathbf{P}_2$ on the coarse mesh $\mathfrak{M}$;

10. &nbsp;&nbsp;&nbsp;&nbsp;Interpolate and get the FE fields $E_{\text{hom}}$, $K_{\text{hom}}$ and $\mathbf{P}_2$ on the fine mesh **m**;

11. &nbsp;&nbsp;&nbsp;&nbsp;Solve macroscopic problem on **m** (Algorithm 3);

12. &nbsp;&nbsp;&nbsp;&nbsp;Compute strain and stress field on **m**;

13. &nbsp;&nbsp;&nbsp;&nbsp;Compute growth rate depending on strain or stress;

14. &nbsp;&nbsp;&nbsp;&nbsp;Compute new growth tensor $\mathbf{F}_g$ on **m**;

15. &nbsp;&nbsp;&nbsp;&nbsp;Project $\mathbf{F}_g$ on the coarse mesh $\mathfrak{M}$;

16. &nbsp;&nbsp;&nbsp;&nbsp;Set $t = t + dt$;

**Algorithm 2: ‘Unit cell’ problems**

**Inputs**: ‘unit cell’ geometry ($l_1, l_2, \theta, w$),
cell wall properties ($E, \nu$),
$P_2$,
growth tensor $\mathbf{F}_g$

**Output**: effective properties $\mathbb{E}_{\text{hom}}$, $K_{\text{hom}}$ and $\mathbf{P}_2$

1. Construct ‘unit cell’ mesh;

2. Solve the ‘unit cell’ problems (67) and (68);

3. Compute the effective properties $\mathbb{E}_{\text{hom}}$ (52), $K_{\text{hom}}$ (53) and $\mathbf{P}_2$ (70).

**Algorithm 3: Macroscopic problem**

**Inputs**: fine mesh **m** of the tissue,
pressure field $P_1$ and growth tensor field $\mathbf{F}_g$ on **m**
homogenized property fields $\mathbb{E}_{\text{hom}}$, $K_{\text{hom}}$ and $\mathbf{P}_2$ on **m**

**Output**: macroscopic displacement field $\mathbf{u}$ on **m**

1. Solve macroscopic problem (64).

For the ‘unit cell’ problems we implement periodic boundary conditions on $\partial Y$ and the corresponding Neumann boundary conditions on $\Gamma$. For the macroscopic problem zero normal displacement and no tangential force are imposed on the lower boundary of the tissue, with zero-force conditions on all other boundaries of $\Omega$. The weak formulations of the ‘unit cell’ problems and the macroscopic problem were implemented using
FreeFEM. The vector components of the displacement fields are represented by P2 finite elements for the ‘unit cell’ problems, while they are P1 elements for the macroscopic problem. We use the sparsesolver solver of FreeFEM to solve both problems.

While the building blocks Algorithms 2 and 3 are implemented in FreeFEM, the orchestrating Algorithm 1 is implemented in Python. In particular, at each time step we solve a relatively high number of ‘unit cell’ problems that are independent of each other, and we use the multiprocessing module of Python to solve them in parallel.

## 5. NUMERICAL SIMULATION RESULTS

In this section we first present the results of the ‘unit cell’ problems, which allow us to compute the effective material properties, and in particular we analyse how the microscopic geometrical parameters as well as growth affect the tissue level properties. Then we validate the multiscale, coupled simulation algorithm on homogeneous isotropic tissues and homogeneous pressure comparing its outcome to the simulation of the microscopic model of the same setup. We also validate the multiscale approach on more general configurations. Namely, we consider gradient fields over the tissue of all the microscopic parameters one by one and, whenever possible, compare the tissue deformation computed using the coupled simulation to the one obtained by simulating the same setup using the microscopic model.

### 5.1. The ‘unit cell’ problems and the macroscopic (homogenized) material properties.
In the case of non-growing cellular solids, previous studies have computed numerically the homogenized elasticity tensor and compared their results to asymptotic expansions and experimental results, see e.g. [34]. We validated our results by quantitatively comparing the dependence of the effective material properties embedded in $E_{hom}$ with the results presented in [34]. Qualitatively, a ‘unit cell’ in the form of a regular hexagon yields an isotropic homogenized medium, see Figure 1d-e. Elongating the unit cell in the $x_2$-direction only reduces the homogenized modulus in the $x_1$-direction. This is in agreement with the measurements on epidermal onion peels realized using a microextensometer setup [33]. The macroscopic (homogenized) Poisson’s ratio $\nu_{12}$ is comparatively large for thin walls and is reduced by elongating the unit cell in the $x_2$-direction, see Figure 1f. The shear modulus turns out to be particularly low for thin walls, and is lower for elongated cells than for regular hexagons, see Figure 1g. All homogenized properties converge to the microscopic properties when the empty space (cell inside) vanishes, i.e. $w/l_2 \to 1$.

In our multiscale analysis we introduced a tensorial material property, $K_{hom}$, that accounts for the contribution of cell pressure to the macroscopic (homogenized) stress tensor. For the present choice of coordinate system (which corresponds to symmetry axes of the unit cell), $K_{hom}$ is diagonal; the two diagonal elements of $K_{hom}$ are shown in Figure 1h-i as functions of relative thickness. A regular hexagonal unit cell yields an isotropic $K_{hom}$, while an elongated hexagon an anisotropic $K_{hom}$, where the anisotropy increases with wall thickness (as it can be seen from the ratio $K_{hom11}/K_{hom22}$).

In case of isotropic growth, neither $E_{hom}$ nor $K_{hom}$ depends on the degree of the growth. In contrast, when growth is anisotropic, the effective material properties change with growth anisotropy. In order to illustrate this, we show in Figure 1j-o the macroscopic (homogenized) material properties as a function of the degree of growth $g$, when there is only growth along axis 1, for a fixed relative wall thickness. The macroscopic (homogenized) Young’s modulus in growth direction increases with $g$, see Figure 1j, while the perpendicular modulus decreases, see Figure 1k. Inversely, the absolute value of $K_{hom}$ in growth direction decreases with $g$, see Figure 1n, while the absolute value of the other element increases, see Figure 1o. The macroscopic (homogenized) Poisson’s ratio $\nu_{12}$ increases with $g$. Finally, the macroscopic (homogenized) shear modulus is non-monotonic.

### 5.2. Validation and sensitivity analysis for tissues with homogeneous material properties.
In order to further validate the multiscale method and the macroscopic model derived from the microscopic description of the growth and elastic deformations, we first consider on the one hand a tissue of $N_x^2/2$, with $N_x = 16$, identical regular hexagonal cells with microscopic geometrical parameters the cell wall length $4\delta l_1$, $4\delta l_2$, the cell wall thickness $4\delta w$, and the angle $\theta$, where $\delta = 1/N_x = 1/16$, $l_1 = l_2 = 1$, $w = 0.05$ and $\theta = \pi/3$, and on the other hand a continuous tissue of the same size. The cell wall material properties, $E = 1$ and $\nu = 0.3$, and the pressure, $P_1 = 0.001$, are homogeneous over the tissue. In Figure 2d-f, we superimposed the initial states of the microscopic cellular structure and of the continuous tissue, and we visually find that the continuous tissue evolved with the macroscopic model corresponds well to the cellular tissue evolved with the microscopic model at $t = 60$, both for the strain- and the stress-based growth law, see Figure 2e-f.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. **Unit cell problem.** a. ‘Unit cell’ examples: isotropic cell shape $l_1/l_2 = 1$ (red box), elongated cell shape $l_1/l_2 = 1.5$ (green box) and $l_1/l_2 = 2$ (blue box) with $w/l_2 = 0.05$ and $\theta = 30$ for all. The color of the box corresponds to the color of lines in other graphs. b. Geometry of the ‘unit cell’ and definition of the geometrical parameters $l_1, l_2, w$ and $\theta$; c. Deformed ‘unit cell’ geometries corresponding to the four elementary solutions of the ‘unit cell’ problems (parameters $E = 1, \nu = 0.35, l_1 = 1, l_2 = 1, w = 0.1, \theta = 30, F_g = I$); **d-i. Dependence of homogenized material properties on the relative wall thickness (no growth).** All moduli are normalized by the cell wall Young modulus $E$ and are shown as a function of relative wall thickness $w/l_2$. **d. Modulus** along axis 1, $E_{hom1}$; **e. Modulus** along axis 2, $E_{hom2}$; **f. First Poisson’s ratio** $\nu_{12}$; **g. Shear modulus**, $G_{12}$; **h-i. Diagonal elements** of the effective material property $K_{hom}$. The off-diagonal elements vanish. **j-o. Dependence of homogenized material properties on the degree of anisotropic growth.** We consider here only growth along axis 1, with corresponding eigenvalue $g$. All moduli are normalized by the cell wall Young modulus $E$ and are shown as a function of $g$. The relative wall thickness is set to $w/l_2 = 0.05$. **j. Modulus** along axis 1, $E_{hom1}$; **k. Modulus** along axis 2, $E_{hom2}$; **l. First Poisson’s ratio** $\nu_{12}$; **m. Shear modulus**, $G_{12}$; **n-o. Diagonal elements** of the effective material property $K_{hom}$. The off-diagonal elements vanish.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIGURE 2. The coupled simulation and its validation. a. The coupled simulation involves in parallel a fine mesh $\mathbf{m}$ and a coarse mesh $\mathcal{M}$ of the tissue; b. Simplified diagram of the coupled simulation algorithm, emphasising the usage of the two meshes; c. One-factor-at-a-time sensitivity of growth rate to model parameters around their reference values ($E = 1, \nu = 0.3, l_1 = l_2 = 1, w = 0.05, \theta = 30, \eta = 1, \tau = 0$), for both stress-based and strain-based growth; d-f. Superimposition of the cellular (cell number, $N_x = 16$) and continuous representations of tissue at $t = 0$ (e) and at $t = 60$ using either strain-based (e) or stress-based (f) growth-law – the smaller shaded area in e-f indicates the initial tissue geometry; g-h. Relative error (difference between coupled simulations and full microscopic model) of the displacement components $u_1$ and $u_2$ as a function of the cell number $N_x$; i-j. Mean growth rate as a function of time, with strain-based (i) or stress-based (j) growth law. The coupled simulation is compared with the microscopic simulation for several values of cell number, $N_x$, and of wall thickness to length ratio, $w/l_2$ – all but $w/l_2$ parameters are the reference parameters.</figcaption></figure>
We varied the number of cells by scaling their size while keeping tissue dimensions unchanged, and considered three values of relative wall thickness, $w/l_2 = 0.05, 0.1$ and $0.2$. We analysed the relative error given by the relative difference of solutions obtained from the coupled macroscopic simulations and microscopic simulations, see (72) in Appendix. Figures 2g-h show the relative error of the displacement components $u_1$ and $u_2$ at the first time step as a function of the number of cells $N_x$. This relative error is higher for thinner walls and it decays with the number of cells. Accordingly, the values of the growth rate are in agreement between microscopic and coupled macroscopic simulations, see Figures 2i-j. The larger error for a tissue of cells with thin walls compared to the error for a tissue of cells with thick walls for the same number of cells relates to the fact that the tissue with thick cell walls has smaller voids and is closer in the approximation to the homogeneous tissue than the tissue with thin cell walls. The mean growth appears constant in time and space, with thinner walls growing faster. This analysis verifies that the coupled simulation scheme enables the efficient computation of the tissue-scale behavior for a given microscopic cellular geometry and cell wall material properties, and provides a good approximation to the microscopic description of the problem.

In order to better understand the effects of microscopic parameters on the macroscopic behavior we performed a one-factor-at-a-time sensitivity analysis around reference values. We assumed the tissue to be homogeneous and we computed the sensitivity of growth rate with respect to each parameter, as defined by equation (73) in Appendix, see Figure 2c. As could be expected, growth rate increases with pressure, $P_1$, and extensibility, $\eta$, while it decreases with yield threshold, $\tau$. Growth rate is less sensitive to yield threshold then to pressure and extensibility. Concerning microscopic geometric parameters (cell size, $l_2$, and wall thickness, $w$) growth is promoted by thinner walls if cell size is constant or by larger cells if wall thickness is constant.

Finally, we note a fundamental difference between the two growth models: while the strain-based growth rate decreases with cell wall Young modulus $E$ and Poisson ratio $\nu$, these microscopic material properties have no effect on the stress-based growth.

**5.3. Validation and predictions based on the dynamics of tissues with heterogeneous material properties.** Now we consider more complex configurations in which material or geometric properties vary spatially, and we use them to further validate the agreement of macroscopic coupled simulations with simulation results for the microscopic model.

Beforehand we perform a parameter scan on larger intervals in the context of homogeneous tissues, and we present in the first column of Figure 3 the growth rate as a function of the parameters for strain- and stress-based growth hypotheses, as predicted by the macroscopic model. The reference situation is highlighted by black dots on each graph. For several parameters the dependence is affine, so that the sensitivity analysis in Figure 2c, which was restricted to small variations around the reference values, already captured the essence of the behavior. We find a non-affine dependence of the growth rate on the cell wall Young modulus for strain based growth, and on cell wall thickness for both growth hypotheses. We also note that the results for the strain-based and stress-based growth rates are identical when the Poisson ratio is equal to zero and all other parameters take their reference values, providing an additional check of the numerical simulation codes.

As next, we consider spatial heterogeneity of each parameter one-by-one, where all other parameters are homogeneous and take the reference value. We consider linear variations of the parameter along the first axis of the tissue, spanning the same interval for which the growth rate variation was presented in the homogeneous context, see the first column of Figure 3. In the second and third columns of Figure 3 we present the superposition of the corresponding tissues at $t = 60$ for the strain-based and stress-based growth hypothesis, respectively. For the five parameters, $E, \nu, P_1, \tau, \eta$, there is a quantitative agreement of tissue shape and size between macroscopic coupled simulations and microscopic simulations. This outcome further validates the multiscale procedure in a context where an additional mechanical stress is induced by spatial differences in growth rate.

Finally, we note that implementing gradients of geometrical properties is not trivial in the microscopic model, whereas such gradients can be easily implemented in the macroscopic coupled simulations. In Figure 3f-g, we present initial and final tissue shape for gradients in thickness, $w$, and cell size, $l_2$, as predicted by the macroscopic two-scale model.

## 6. CONCLUSIONS

We formulated a microscopic model for the growth of plant tissues and derived the corresponding macroscopic equations, when the limit of the ratio between cell size and tissue size $\delta$ tends to zero. Both methods, the formal asymptotic expansion and rigorous two-scale convergence, are used to revise the macroscopic equations. Passing to the limit in the nonlinear equations, resulted from the multiplicative decomposition of the deformation gradient into elastic and growth parts, requires the strong convergence of $\{F_{g,\delta}\}$, the proof of which was
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Validation and illustration with gradients of material/geometric properties.
First column: One factor-at-a-time parameter scan around the reference model for the homogeneous macroscopic coupled simulation. The results for the reference parameter values are indicated by black dots. Second and third columns: Superimposition of the deformed reference tissues obtained from microscopic (excluding f and g) and from macroscopic coupled simulations, where all parameters are homogeneous, except for the parameter of focus in the row, which varies linearly with first coordinate (in initial geometry) between the minimum and maximum value of the x-axis in first column. The shadowed area behind indicates the initial geometry of the tissue. In all simulations $t_{max} = 60$, $dt = 1$; the number of cells is $N_x = 16$ for microscopic simulations. The second and third columns correspond to strain-based and stress-based growth hypotheses, respectively. The parameters varied in each row are: a. Cell-wall Young modulus, $E$; b. Cell-wall Poisson ratio, $\nu$; c. Pressure, $P_1$; d. Strain/stress threshold value, $\tau$; e. Extensibility, $\eta$; f. Wall thickness, $w$; g. Cell size, $l_1 = l_2$.</figcaption></figure>
the main technical step in the analysis. To show the strong convergence for the sequence of solutions of the microscopic problem and uniqueness for the macroscopic problem we used assumptions on the boundedness of the deformation gradient, which may be shown using the regularity results for solutions of the linear elasticity equations. The multiscale analysis performed here will also apply to other multiscale models for the growth of biological tissues.

The macroscopic two-scale model comprises the coupled system of equations of linear elasticity for the displacement and nonlinear ordinary differential equations for the growth tensor, that are additionally coupled to the ‘unit cell’ problems describing cell-scale properties and behaviour. This approach notably enables mapping the parameters of microscopic models that account for cell level characteristics to those of the macroscopic models defined on the tissue level, helping to bridge two relatively separate worlds in current models of plant growth. We implemented both models (microscopic and macroscopic) numerically and compared them quantitatively. The macroscopic coupled model requires less computational time for its numerical solution and can be readily adapted to describe spatial variations of parameters. Future work will address more complex spatial patterns of properties, such as accounting for variations in material properties across the cell wall, or additional physics, such as feedback from stress on material properties or hydraulics of water movement during growth. More generally, homogenization and the corresponding multiscale models appear as promising approaches to quantitatively describe morphogenesis.

## APPENDIX A. COMPONENTS OF THE ELASTIC TENSOR

The elasticity tensor in two dimensions is implemented as a six-component structure, where every component varies spatially,

$$ 
\begin{aligned}
\mathbb{E} &= \{E_0, E_1, E_2, E_3, E_4, E_5\} \\
&= \{E_{1111}, E_{2222}, E_{1122}, E_{1212}, E_{1112}, E_{2212}\}.
\end{aligned}
$$

For an isotropic material in 2 dimensions, or in 3 dimensions with plane stress conditions, the components of the elasticity tensor using Lamé's parameters $\lambda$ and $\mu$ are given by

$$ \mathbb{E} = \{2\mu + \lambda, 2\mu + \lambda, \lambda, \mu, 0, 0\}, \quad \text{with} \quad \mu = \frac{E}{2(1 + \nu)}, \quad \lambda = \frac{E\nu}{1 - \nu^2}, $$

where $E$ and $\nu$ are Young's modulus and Poisson's ratio respectively.

Given an elasticity tensor $\mathbb{E}$ in 2 spatial dimensions, for an orthotropic material with symmetry axes along the two Cartesian axes ($E_{1112} = E_{2221} = 0$), its material properties, such as the Young moduli $E_1$ and $E_2$ in the symmetry directions as well as the Poisson ratio $\nu_{12}$ and shear modulus $G_{12}$, can be computed as

$$ E_1 = E_{1111} - \frac{E_{1122}^2}{E_{2222}}, \quad E_2 = E_{2222} - \frac{E_{1122}^2}{E_{1111}}, \quad \nu_{12} = \frac{E_{1122}}{E_{2222}}, \quad G_{12} = E_{1212}. $$

## APPENDIX B. ERROR AND SENSITIVITY DEFINITIONS

For the quantity $a^\delta$ defined in the microscopic domain (tissue) $\Omega^\delta$ and its corresponding macroscopic quantity $A$ defined in the macroscopic domain (tissue) $\Omega$, we define the relative error $e^\delta$ as

$$ (72) \qquad e^\delta(a^\delta, A) = \sqrt{\frac{\int_{\Omega^\delta} (a^\delta - A)^2 dx}{\int_{\Omega^\delta} A^2 dx}}. $$

The sensitivity $\phi(X, x)$ around the reference model of an output quantity $X$ to the value of parameter $x$ is defined as

$$ (73) \qquad \phi(X, x) = \frac{x_0}{X(x_0)} \frac{\partial X}{\partial x}\bigg|_{x=x_0}, $$

where $x_0$ is the reference value of the input parameter $x$ and $X(x_0)$ is the output value in the reference model (all input parameters are equal to their reference value). The sensitivities $\phi$ are normalized, non-dimensional quantities, and thus sensitivities of the output to different input parameters can be compared.

## APPENDIX C. TWO-SCALE CONVERGENCE AND PERIODIC UNFOLDING OPERATOR

We recall the definition and some properties of the two-scale convergence and periodic unfolding operator.

**Definition C.1** (Two-scale convergence). [3, 41] A sequence $\{u^\delta\}$ in $L^p(\Omega)$, with $1 < p < \infty$, is two-scale convergent to $u \in L^p(\Omega \times Y)$ if for any $\phi \in L^q(\Omega; C_{\text{per}}(Y))$, with $1/p + 1/q = 1$,

$$ \lim_{\delta \to 0} \int_{\Omega} u^{\delta}(x) \phi(x, x/\delta) dx = \frac{1}{|Y|} \int_{\Omega \times Y} u(x, y) \phi(x, y) dy dx. $$
**Theorem C.2.** [5, 40] Let {$v^\delta$} $\subset L^2(\Gamma^\delta)$ satisfy $\delta\|v^\delta\|_{L^2(\Gamma^\delta)}^2 \le C$, then there exists a two-scale limit $v \in L^2(\Omega; L^2(\Gamma))$ such that, up to a subsequence, $v^\delta$ two-scale converges to $v \in L^2(\Omega \times \Gamma)$ in the sense that

$$ \lim_{\delta \to 0} \delta \int_{\Gamma^\delta} v^\delta(x) \phi(x, x/\delta) d\gamma^\delta = \frac{1}{|Y|} \int_{\Omega \times \Gamma} v(x, y) \phi(x, y) d\gamma_y dx, \quad \text{for any } \phi \in C_0(\Omega; C_{\text{per}}(Y)). $$

**Lemma C.3.** (i) If {$u^\delta$} is bounded in $L^2(\Omega)$, there exists a subsequence (not relabelled) such that $u^\delta \to u$ two-scale as $\delta \to 0$ for some function $u \in L^2(\Omega \times Y)$.

(ii) If $u^\delta \to u$ weakly in $H^1(\Omega)$ then $u^\delta \to u$ and $\nabla u^\delta \to \nabla u + \nabla_y u_1$ two-scale, where $u_1 \in L^2(\Omega; H_{\text{per}}^1(Y)/\mathbb{R})$.

To define the periodic unfolding operator, let $[z]_Y$ for any $z \in \mathbb{R}^d$ denote the unique combination $\sum_{i=1}^d k_i b_i$ with $k \in \mathbb{Z}^d$, such that $z - [z]_Y \in Y$, see e.g. [11, 13]

**Definition C.4.** Let $p \in [1, \infty]$ and $\phi \in L^p(\Omega)$. The unfolding operator $T^\delta$ is defined by

$$ T^\delta(\phi)(x, y) = \begin{cases} \phi(\delta[x/\delta]_Y + \delta y) & \text{for a.e. } (x, y) \in (\Omega \setminus \Lambda_\delta) \times Y, \\ 0 & \text{for a.e. } x \in \Lambda_\delta, y \in Y, \end{cases} \quad \text{and} \quad T^\delta(\phi) \in L^p(\Omega \times Y). $$

For $\psi \in L^p(\Gamma^\delta)$ the boundary unfolding operator $T^\delta$ is defined by

$$ T^\delta(\psi)(x, y) = \begin{cases} \psi(\delta[x/\delta]_Y + \delta y) & \text{for a.e. } (x, y) \in (\Omega \setminus \Lambda_\delta) \times \Gamma, \\ 0 & \text{for a.e. } x \in \Lambda_\delta, y \in \Gamma, \end{cases} \quad \text{and} \quad T^\delta(\psi) \in L^p(\Omega \times \Gamma). $$

For $\psi \in L^p(\Omega^\delta)$ the unfolding operator $T^\delta$ is defined in (24).

Notice that in the main text we use the same notation $T^\delta$ for all three types of unfolding operator.

**Proposition C.5** ([12]). Let {$\psi^\delta$} be a bounded sequence in $L^p(\Omega)$ for some $1 < p < \infty$. Then the following assertions are equivalent:

(i) {$T^\delta(\psi^\delta)$} converges weakly to $\psi$ in $L^p(\Omega \times Y)$; (ii) {$\psi^\delta$} converges two-scale to $\psi$, $\psi \in L^p(\Omega \times Y)$.

We have the following properties of the periodic unfolding operator and the boundary unfolding operator:

$$ \begin{align} T^\delta(F(u, v)) &= F(T^\delta(u), T^\delta(v)), \quad T^\delta(v(t, x/\delta)) = v(t, y), \quad x \in \Omega^\delta, y \in Y_w \text{ or } y \in \Gamma, t \in (0, T), \tag{74} \\ \langle T^\delta(v), T^\delta(u) \rangle_{\Omega_T \times Y_w} &= |Y| \langle v, u \rangle_{\Omega_T^\delta} - |Y| \langle v, u \rangle_{\Lambda_{\delta,T}}, \quad \|T^\delta(\phi)\|_{L^p(\Omega_T \times Y_w)} \le |Y|^{\frac{1}{p}} \|\phi\|_{L^p(\Omega_T^\delta)}, \notag \\ \|T^\delta(\psi)\|_{L^p(\Omega_T \times \Gamma)} &\le \delta^{\frac{1}{p}} |Y|^{\frac{1}{p}} \|\psi\|_{L^p(\Gamma_T^\delta)} \le C (\|\psi\|_{L^p(\Omega_T^\delta)} + \delta \|\nabla \psi\|_{L^p(\Omega_T^\delta)}), \notag \end{align} $$

for $u, v \in L^2((0, T) \times \Omega^\delta)$ or $u, v \in L^2((0, T) \times \Gamma^\delta)$, $\phi \in L^p((0, T) \times \Omega^\delta)$, $\psi \in L^p(0, T; W^{1,p}(\Omega^\delta))$ and $F$ is any linear or nonlinear function, see e.g. [12, 11, 13].

**Lemma C.6** ([12]). (i) If $\phi \in L^p(\Omega)$, then $T^\delta(\phi) \to \phi$ strongly in $L^p(\Omega \times Y)$, for $1 \le p < \infty$.

(ii) Let {$\psi^\delta$} $\subset L^p(\Omega)$, with $\psi^\delta \to \psi$ strongly in $L^p(\Omega)$, then $T^\delta(\psi^\delta) \to \psi$ strongly in $L^p(\Omega \times Y)$.

## NUMERICAL SIMULATION CODES

Numerical simulation codes for both microscopic model (8), (11) and macroscopic model (54), (55) can be found under https://gitlab.inria.fr/akiss1/planthom.

## ACKNOWLEDGMENTS

We would like to thank the Newton Institute of Mathematical Sciences (INI) and International Centre for Mathematical Sciences (ICMS) for organising the research programme 'Growth, form and self-organisation' and the workshop 'Growth, form and self-organisation in living systems' at which this work was started. We also acknowledge the support of the Centre Blaise Pascal's IT test platform operated with SIDUS [48] at ENS de Lyon for the possibility to test our code on the center's computers.
REFERENCES

[1] S. AGMON, A. DOUGLIS, AND L. NIRENBERG, *Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions*, i, Comm. Pure Appl. Math., 12 (1959), pp. 623-727.

[2] S. AGMON, A. DOUGLIS, AND L. NIRENBERG, *Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions*, ii, Comm. Pure Appl. Math., 17 (1964), pp. 35-92.

[3] G. ALLAIRE, *Homogenization and two-scale convergence*, SIAM Journal of Mathematics and Analysis, 23 (1992), pp. 1482-1518.

[4] G. ALLAIRE, *Shape Optimization by the Homogenization Method*, Springer-Verlag New York, 2002.

[5] G. ALLAIRE, A. DAMLAMIAN, AND U. HORNUNG, *Two-scale convergence on periodic surfaces and applications*, in Proceedings of the International Conference on Mathematical Modelling of Flow Through Porous Media, A. Bourgeat, C. Carasso, and A. Mikelic, eds., World Scientific, Singapore, (1996), pp. 15-25.

[6] H. ALT, *Linear Functional Analysis. An Application-Oriented Introduction*, Springer, 2016.

[7] B. BOZORG, P. KRUPINSKI, AND H. JÖNSSON, *Stress and strain provide positional and directional cues in development*, PLOS Computational Biology, 10 (2014), p. e1003410.

[8] B. BOZORG, P. KRUPINSKI, AND H. JÖNSSON, *A continuous growth model for plant tissue*, Phys. Biol., 13 (2016), p. 065002.

[9] A. CIANCHI AND V. MAZYA, *Global boundedness of the gradient for a class of nonlinear elliptic systems*, Arch. Rational Mech. Anal., 212 (2014), pp. 129-177.

[10] P. CIARLET, *Mathematical elasticity. Volume I: Three-dimensional elasticity*, North-Holland, 1988.

[11] D. CIORANESCU, A. DAMLAMIAN, P. DONATO, G. GRISO, AND R. ZAKI, *The periodic unfolding method in domains with holes*, SIAM Journal of Mathematics and Analysis, 44 (2012), pp. 718-760.

[12] D. CIORANESCU, A. DAMLAMIAN, AND G. GRISO, *The periodic unfolding method in homogenization*, SIAM Journal of Mathematics and Analysis, 40 (2008), pp. 1585-1620.

[13] D. CIORANESCU, A. DAMLAMIAN, AND G. GRISO, *The Periodic Unfolding Method*, Springer, 2019.

[14] D. CIORANESCU AND P. DONATO, *An Introduction to Homogenization*, Oxford University Press, 1999.

[15] F. CORSON, M. ADDA-BEDIA, AND A. BOUDAOUD, *In silico leaf venation networks: growth and reorganization driven by mechanical forces*, Journal of Theoretical Biology, 259 (2009), pp. 440-448.

[16] E. DAVOLI, C. GAVIOLI, AND V. PAGLIARI, *A homogenization result in finite plasticity*, arXiv::2204.09084v2, (2022).

[17] E. DAVOLI, K. NIK, AND U. STEFANELLI, *Existence results for a morphoelastic model*, J Appl Math Mech, (2022).

[18] G. DUVAUT AND J. LIONS, *Inequalities in Mechanics and Physics*, Springer, 1976.

[19] R. DYSON, L. BAND, AND O. JENSEN, *A model of crosslink kinetics in the expanding plant cell wall: Yield stress and enzyme action*, Journal of Theoretical Biology, 307C (2012), pp. 125-136.

[20] R. DYSON AND O. JENSEN, *A fibre-reinforced fluid model of anisotropic plant cell growth*, Journal of Fluid Mechanics, 655 (2010), pp. 472-503.

[21] L. EVANS, *Partial Differential Equations*, AMS, 2010.

[22] A. GORIELY, *The Mathematics and Mechanics of Biological Growth*, Springer, 2017.

[23] A. GORIELY AND M. B. AMAR, *On the definition and modeling of incremental, cumulative, and continuous growth laws in morphoelasticity*, Biomechanics and Modeling in Mechanobiology, 6 (2007), pp. 289-296.

[24] O. HAMANT, M. HEISLER, H. JÖNSSON, P. KRUPINSKI, M. UYTTEWAAL, P. BOKOV, F. CORSON, P. SAHLIN, A. BOUDAOUD, E. MEYEROWITZ, Y. COUDER, AND J. TRAAS, *Developmental patterning by mechanical signals in arabidopsis*, Science, 322 (2008), pp. 1650-1655.

[25] F. HECHT, *New development in FreeFem++*, J. Numer. Math., 20 (2012), pp. 251-265.

[26] N. HERVIEUX, M. DUMOND, A. SAPALA, A. ROUTIER-KIERZKOWSKA, D. KIERZKOWSKI, A. ROEDER, R. SMITH, A. BOUDAOUD, AND O. HAMANT, *A mechanical feedback restricts sepal growth and shape in arabidopsis*, Current Biology, 26 (2016), pp. 1019-1028.

[27] U. HORNUNG AND W. JÄGER, *Diffusion, convection, adsorption and reaction of chemicals in porous media*, J. Differential Equations, 92 (1991), pp. 199-225.

[28] R. HUANG, A. BECKER, AND I. JONES, *A finite strain fibre-reinforced viscoelasto-viscoplastic model of plant cell wall growth*, Journal of Engineering Mathematics, 95 (2015), pp. 121-154.

[29] V. V. JIKOV, S. M. KOZLOV, AND O. A. OLEINIK, *Homogenization of Differential Operators and Integral Functionals*, Springer-Verlag, 1994.

[30] C. KENIG, F. LIN, AND Z. SHEN, *Homogenization of elliptic systems with Neumann boundary conditions*, J Americal Math Society, 26 (2013), pp. 901-937.

[31] J. KHADKA, J.-D. JULIEN, AND K. ALIM, *Feedback from tissue mechanics self-organizes efficient outgrowth of plant organ*, Biophysical J, 117 (2019), pp. 1995-2004.

[32] J. LOCKHART, *An analysis of irreversible plant cell elongation*, Journal of Theoretical Biology, 8 (1965), pp. 264-275.

[33] M. MAJDA, N. TROZZI, G. MOSCA, AND R. SMITH, *How Cell Geometry and Cellular Patterning Influence Tissue Stiffness*, International Journal of Molecular Sciences, 23 (2022), p. 5651.

[34] S. MALEK AND L. GIBSON, *Effective elastic properties of periodic hexagonal honeycombs*, Mechanics of Materials, 91 (2015), pp. 226-240.

[35] R. M. H. MERKS, M. GURAVAGE, D. INZE, AND G. T. S. BEEMSTER, *VirtualLeaf: an open-source framework for cell-based modeling of plant tissue growth and development*, Plant Physiology, 155 (2011-02), pp. 656-666.

[36] A. MIELKE AND U. STEFANELLI, *Linearized plasticity is the evolutionary $\gamma$-limit of finite plasticity*, J. Eur. Math. Soc., 15 (2013), pp. 923-948.

[37] J. NEČAS, *Direct Methods in the Theory of Elliptic Equations*, Springer, 2012.

[38] P. NEFF, *On Korn's first inequality with non-constant coefficients*, Proc Royal Society of Edinburgh, 132A (2002), pp. 221-243.
[39] P. NEFF, *Local existence and uniqueness for quasistatic finite plasticity with grain boundary relaxation*, Quarterly of Applied Mathematics, 63 (2005), pp. 88–116.

[40] M. NEUSS-RADU, *Some extensions of two-scale convergence*, C. R. Math. Acad. Sci. Paris, 332 (1996), pp. 899–904.

[41] G. NGUETSENG, *A general convergence result for a functional related to the theory of homogenization*, SIAM Journal of Mathematics and Analysis, 20 (1989), pp. 608–623.

[42] O. A. OLEINIK, A. S. SHOMAEV, AND G. A. YOSIFIAN, *Mathematical Problems in Elasticity and Homogenization*, North-Holland, 1992.

[43] J. A. OLIVEIRA, J. PINHO-DA CRUZ, AND F. TEIXEIRA-DIAS, *Asymptotic homogenisation in linear elasticity. Part II: Finite element procedures and multiscale applications*, Computational Materials Science, 45 (2009), pp. 1081–1096.

[44] H. OLIVERI, J. TRAAS, C. GODIN, AND O. ALI, *Regulation of plant cell wall stiffness by mechanical stress: a mesoscale physical model*, Journal of Mathematical Biology, 78 (2019), pp. 625–653.

[45] J. PINHO-DA CRUZ, J. A. OLIVEIRA, AND F. TEIXEIRA-DIAS, *Asymptotic homogenisation in linear elasticity. Part I: Mathematical formulation and finite element modelling*, Computational Materials Science, 45 (2009), pp. 1073–1080.

[46] W. POMPE, *Korn's first inequality with variable coefficients and its generalization*, Comment. Math. Univ. Carolin, 44 (2003), pp. 57–70.

[47] M. PTASHNYK AND B. SEGUIN, *The impact of microfibril orientations on the biomechanics of plant cell walls and tissues: modelling and simulations*, Bulletin of Mathematical Biology, 78 (2016), pp. 2135–2164.

[48] E. QUEMENER AND M. CORVELLEC, *SIDUS—the Solution for Extreme Deduplication of an Operating System*, The LINUX Journal, January (2014).

[49] E. RODRIGUEZ, A. HOGER, AND A. MCCULLOCH, *Stress-dependent finite growth in soft elastic tissue*, J Biomechanics, 27 (1994), pp. 455–467.

[50] E. SANCHEZ-PALENCIA, *Non-homogeneous media and vibration theory*, Springer-Verlag Berlin, 1980.

[51] Z. SHEN, *Boundary estimates in elliptic homogenization*, Anal. PDE, 3 (2017), pp. 653–694.

[52] J. SIMON, *Compact sets in the space $l^p(0, t; b)$*, Annali di Matematica Pura ed Applicata volume, 146 (1986), pp. 65–96.

[53] E. SMITHERS, J. LUO, AND R. DYSON, *Mathematical principles and models of plant growth mechanics: from cell wall dynamics to tissue morphogenesis*, Journal Of Experimental Botany, 70 (2019), pp. 3587–3600.

[54] R. VANDIVER AND A. GORIELY, *Tissue tension and axial growth of cylindrical structures in plants and elastic tissues*, EPL (Europhysics Letters), 84 (2008), p. 58004.