Geometry of adipocyte packing in subcutaneous tissue contributes to nonlinear tissue properties captured through a Gaussian process surrogate model

Abstract

Subcutaneous tissue mechanical response is governed by the geometry and mechanical properties at the microscale and drives physiological and clinical processes such as drug delivery. Even though adipocyte packing is known to change with age, disease, and from one individual to another, the link between the geometry of the packing and the overall mechanical response of adipose tissue remains poorly understood. Here we create 1200 periodic representative volume elements (RVEs) that sample the possible space of Laguerre packings describing adipose tissue. RVE mechanics are modeled under tri-axial loading. Equilibrium configuration of RVEs is solved by minimizing an energetic potential that includes volume change contributions from adipocyte expansion, and area change contributions from collagen foam stretching. The resulting mechanical response across all RVE samples is interpolated with the aid of a Gaussian process (GP), revealing how the microscale geometry dictates the overall RVE mechanics. For example, increase in adipocyte size and increase in sphericity lead to adipose tissue softening. We showcase the use of the homogenized model in finite element simulations of drug injection by implementing a Blatz–Ko model, informed by the GP, as a custom material in the popular open-source package FEBio. These simulations show how microscale geometry can lead to vastly different injection dynamics even if the constituent parameters are held constant, highlighting the importance of characterizing individual's adipose tissue structure in the development of personalized therapies.

---

1 Introduction

Subcutaneous tissue, underneath the epidermis and dermis layers of the skin, is characterized by the presence of adipocyte cells engulfed in a collagen foam.¹ Subcutaneous tissue mechanics play a key role in many physiological and clinical settings such as drug delivery,² especially due to the rising adoption of auto-injector devices.³ Similar to other biological tissues, the subcutaneous tissue mechanics are governed by the geometry and mechanical properties at the microscale. In the case of subcutaneous tissue, the overall mechanical properties stem from the individual properties of adipocytes, collagen, and, to a lesser extent, other constituents like GAGs, elastin, among others.^1,4,5 Individual constituent properties are then nonlinearly coupled through the geometry of the adipocyte packing. At the macroscopic scale, on the order of centimeters, the tissue can be described within a continuum mechanics framework using the poroelastic theory.^6,7 While there have been experimental studies at the macroscale to characterize adipose tissue mechanics,⁵ imaging an structure studies to get the geometry of adipocytes,⁸ and characterization of individual cell mechanics,⁹ little has been done to connect the two scales.^1,10 This is particularly important because in many cases it is impossible to get enough tissue to do mechanical testing, but it is feasible to image the structure in small biopsies.¹¹ Additionally, there are changes in microstructure of adipose tissue with age, sex, and disease.^8,12,13 Thus, structure–function models relating adipocyte packing to macroscale mechanical behavior can shed light on adipose tissue disease and improve design of drug delivery devices.

At the macroscale, subcutaneous tissue undergoes large deformations, including large volume changes, especially during drug delivery, and has to be modeled as a poro-hyperelastic material.^2,14 For example, experimental work on subcutaneous tissues from pigs and humans has shown that adipose tissue can deform up to ∼200% before failing under uniaxial tension.^15,16 The response is nonlinear, consistent with the presence of a collagen foam which shows strain stiffening.^1,16 The mechanical behavior of adipose tissue has been shown to be isotropic, consistent with an isotropic packing of adipocytes at the microscale.^17,18 During drug injection, volumetric strains can achieve up to ∼50%. Energy dissipation in adipose tissue can be captured within the poro-hyperelastic framework through the fluid flow in the interstitial space.^19,20 Additional energy dissipation mechanisms such as damage¹⁶ and viscoelasticity^4,5,21 can become relevant under certain loading conditions but are not considered in the present study.

Adipocyte packing geometry is a key indicator of age, sex, and disease which, among other effects, is reflected in changes of the tissue biomechanics.^9,22 Adipocyte volume and number increase concomitantly with age until adulthood, when these metrics stabilize.^23,24 However, a hallmark of metabolic disease is adipocyte hypertrophy, i.e. increased adipocyte size, as opposed to hyperplasia (increase in adipocyte number while keeping individual cell volume constant).⁸ Adipocyte physiology is strongly regulated through sex hormones, which leads to notable differences between males and females with regards to adipose tissue.^24,25 For example, pre-menopausal women tend to show increased subcutaneous adipose tissue storage compared to visceral adipose tissue, whereas the opposite is true in men. Interestingly, subcutaneous tissue is less prone to contribute to insulin resistance and metabolic disease compared to visceral adipose tissue.²⁶ Adipocyte sphericity has been quantified to a much lesser extent than adipocyte size, but it is hypothesized that dysfunctional adipose tissue remodeling can lead to a loss of sphericity of adipocytes.¹³

This study presents a multiscale model for subcutaneous tissue that can help us understand its homogenized mechanical behavior based on its microstructure. The microstructure of subcutaneous tissue consists of a Laguerre packing of adipocytes^17,27 in a collagen foam.¹ To derive the homogenized behavior of subcutaneous tissue, multiple representative volume elements (RVE) of the subcutaneous tissue microstructure were generated by sampling the geometry space of Laguerre packings.^28,29 Using a Gaussian process (GP) model,³⁰ we interpolated the response across RVEs to understand the effect of geometry on the homogenized mechanical response. We showcase the utility of the homogenized model by conducting finite element simulations on the open-source platform FEBio,³¹ using a custom implementation of the Blatz–Ko material model,³² informed by the GP. The multiscale model developed in this study has several potential applications. For instance, predicting the mechanical behavior of subcutaneous tissue during drug injection can be used to optimize auto-injector devices for specific patient populations.³³

2 Methods

2.1 Microstructural model and RVE generation

Subcutaneous tissue is composed of adipocytes in a collagen foam. Cellular packings have often been described with Voronoi tessellations.^34,35 However, Voronoi packings are too restrictive to capture the geometry of cellular tissue such as adipose tissue.³⁶ Laguerre packings are more general, they encompass Voronoi tessellations but, moreover, are able to capture arbitrary packings of convex polyhedra.³⁷ The microstructure of a Laguerre tessellation can be statistically described assuming the cells have log-normal distributions³⁸ for the equivalent diameter and sphericity.³⁷ The mean equivalent diameter³⁹mean sphericity,¹³and the standard deviation of both (assuming a log-normal distribution), are the four parameters that define the geometry of the RVE.⁴⁰ In eqn (1), (2)V^(c) and A^(c) are the volume and surface area of a cell respectively. The mean equivalent diameter 〈d〉 and mean sphericity 〈S〉 will simply be denoted as d and S respectively unless otherwise specified.

The operation 〈•〉 in eqn (1) and (2) denotes averaging and is defined as

where n is the number of cells in a given RVE.

Adipocyte equivalent diameters have been observed to range from 34 μm to 126 μm^8,18,41,42 based on histology images. The size of the adipocytes has been observed to have a strong correlation with BMI and fat content.^22,43 The value for sphericity has been observed to range from 0.78 to 0.85,¹³ although there is limited information on values for this parameter.

Neper,⁴⁴ an open-source software package for polycrystal generation allows us to create RVEs based on the four parameters that define the geometry. Based on the Latin hypercube sampling algorithm, 400 combinations of d, d_std, S and S_std were generated, where (•)_std stands for the standard deviation of the given geometric variable. The boundaries for these parameters, shown in Table 1, were chosen such that they will allow us to explore within the values reported in the literature but we extend the limits to unreported results to capture a much wider range of geometries. Based on these parameters, fully periodic RVEs with side lengths of l₀ = 0.5 mm were created. This size allows to capture the behavior of subcutaneous tissue, as smaller RVEs would not be statistically representative. In the iterative creation of RVEs in Neper, discrepancies may arise between the input geometry parameters and the final parameters of the RVE. For the results below, the d, d_std, S, and S_std are recalculated using eqn (1) and (2) after the generation of each RVE.

Table 1 Range for geometry parameters determining the Laguerre packings

Parameter	Min. value	Max. value
d [μm]	25	140
d std [μm]	1	20
S	0.75	0.95
S std	0.001	0.04

---

In addition, to consider cell allocation and the significant impact of spatial arrangement on the nonlinearity of the structural behavior of loaded cells, the seeding point of the RVEs were varied three times while keeping the four geometric inputs constant. This approach takes into account research findings that highlight the pronounced influence of spatial arrangement on the constitutive behavior of individual cellular components.^9,45 As a result, a total of 1200 RVEs were generated to capture the range of possible configurations. After this, all the RVEs with ρ greater than 32 000 cells per mm³ were discarded as they are computationally expensive to run and not physiologically relevant. Based on^8,18,41,42 we should expect the cell density to vary from 1000 to 8500 cells per mm³ for subcutaneous tissue. Nevertheless, we extend the study to account for RVEs with higher ρ, as mentioned, up to 32 000 cells per mm³. Fig. 1 shows possible configurations of the microstructure of subcutaneous tissue by varying a single parameter at a time with respect to their mean value.

Fig. 1 RVEs based on Table 1 by setting all parameters to their mean (〈d〉, 〈d_std〉, 〈S〉, 〈S_std〉) except: (a) d = 1.2 〈d〉, (b) d_std = 1.2 〈d_std〉, (c) S = 0.99, (d) S_std = 1.2 〈S_std〉, (e) d = 0.8 〈d〉, (f) d_std = 0.8 〈d_std〉, (g) S = 0.8 〈S〉, (h) S_std = 0.8 〈S_std〉.

2.2 Elastic behavior

The changes in cell volume are associated with adipocyte deformation, whereas the deformation of faces of the polyhedra across cells can be associated with the deformation of the collagen foam.⁴⁶ Thus, the energy associated with the RVE can be decomposed into two energy components

ψ = ψ^v + ψ^a, (4)

where ψ^v is the energy associated with the changes in the volume of the adipocytes given as

ψ^v = k_v〈(J^(c) − 1)²〉, (5)

and ψ^a is associated with the changes in the area of the collagen sheets enclosing the adipocytes

ψ^a = k_a〈(J^(c)_a − 1)²〉. (6)

J^(c) = V_x/V_X and J^(c)_a = A_x/A_X are the relative volume change and area change respectively of each cell c where the subscript X is the cell property in the reference configuration and the subscript x is the cell property at the deformed configuration.

The material elastic properties, k_v and k_a, capture the response of the material to changes in volume and area, respectively. Since the effect of geometry is the focus of this study, these properties are kept constant. The chosen values for k_v and k_a are 0.3 kPa and 0.7 kPa, respectively, which produce realistic stress values for subcutaneous tissue under a hydrostatic pressure load.¹⁴ However, further research is needed to determine more appropriate values for these material properties.

The Cauchy stress can be calculated in terms of its principal components by taking derivatives with respect to the principal stretches of RVE deformation,⁴⁷

where λ_i are the principal stretches of the deformation and J = λ₁λ₂λ₃ is the overall volume change. We can apply the central difference approximation such that eqn (7) can be approximated astaking ε = 10⁻⁶ in this study. The reason for the numerical approximation of the stress is that we can only compute the energies through evaluation of a forward equilibrium simulation and thus the derivatives (7) cannot be directly evaluated. Nonetheless, we can run simulations for different λ_i and compute the energy, then estimate the stress numerically. The notation σ^e_i is chosen to link this stress to a poroelastic framework.² In poroelasticity, the total stress σ = σ^e − pI, where p is the fluid pressure and σ^e is the elastic stress, associated with the deformation of the material only.

2.3 Simulation of RVE deformation

In order to study the volume change of subcutaneous tissue under drug delivery injections, equi-triaxial periodic boundary conditions (PBCs) were applied to the RVEs. This is an important consideration from an application standpoint as it mimics hydrostatic pressure loading which occurs at the time fluid is injected in the subcutaneous tissue space. Such a deformation has the deformation gradientwhich was applied to all RVEs. Stretches were gradually increased λ = 1.0, 1.05, 1.10, 1.15, and 1.20 and equilibrium reached at each deformation before moving to the next level of deformation. We consider a poroelastic framework.² Assuming zero mixture normal traction then the observed stretches can be seed as the result of uniform pressure σ^e = pI.

With the consideration of tri-axial PBC and isotropy of the material we can simplify the notation to σ^e = σ^e₁ = σ^e₂ = σ^e₃ and λ = λ₁ = λ₂ = λ₃ and drop the i index.

Dynamic relaxation was employed to obtain the equilibrium configurations.⁴⁸ To ensure PBCs while otherwise allowing nodes on the boundary to move freely, minimization of (4) was done by iteratively modifying nodal positions with a force stemming from the variation of the energy, plus a constraint force,

In eqn (10), xⁱ are the current nodal coordinates in a matrix of size M × 3, with M the number of nodes; dt is a pseudo-time to reach the equilibrium solution iteratively set empirically to 0.0015, and F^PBC are the forces needed to impose periodic boundary conditions

such that they are equal to 0 at equilibrium. n^p and n^q are coupled sets of periodic nodes of the RVE that define the coordinates of the end-points of a boundary edge where p and q are index numbers that can range from 1 to the number of nodes at the boundary. For instance, n^p returns the global nodal index and n^p_j returns the coordinates of said node. Θ^pq_j is the periodicity vector associated with the node pair n^p and n^q specifying in which direction the nodes must satisfy the PBC constraint. The components of Θ^pq_j can only be −1, 0 or 1. This information is extracted from the Neper file generated when creating the RVEs. β^PBC is a penalty factor. We set β^PBC = 100 after trial and error. The higher the β^PBC the better the PBC constraint is enforced, but the slower convergence becomes. We found empirically that β^PBC = 100 satisfies the PBC and has a convergence that we could tolerate with our computational budget. For a particular stretch, the change in length of the RVE is calculated as dl = l₀(λ − 1).

To calculate the partial derivatives ∂ψ^v/∂x and ∂ψ^a/∂x, we employ the convex hull algorithm⁴⁹ available in ref. 50 to determine the volume and area of each cell based on the nodes that compose them. The dynamic relaxation eqn (10) is then iteratively solved until the norm of dt(∂ψ^v/∂x + ∂ψ^a/∂x + F^PBC) reaches a set tolerance value (10⁻⁵). As shown in Fig. 2, the volume and area of each cell changes after deformation. While finding the final nodal position, we allow the convexity of each cell to drop, if necessary, with respect to the original Laguerre packing, for which all cells are convex polyhedra.

Fig. 2 (a) Cell volume plotted over the initial and deformed RVE, (b) corresponding distribution of relative volume change of the cells J^(c), (c) cell area in the reference and deformed RVEs, and (d) distribution of relative area change J^(c)_a over the cells after a deformation of λ = 1.15 for a random RVE characterized by d = 91.8 μm, d_std = 12.9 μm, S = 0.89, S_std = 0.015 with n = 291.

2.4 Gaussian process interpolation

After filtering the RVEs with acceptable density ρ, data was generated with 45% of the data points (N = 2359) after randomly mixing the data set. The remaining 55% of the data points were set aside for model validation, resulting in an RMSE value of 0.4556 kPa. To better understand the relationship between geometry and mechanical response, a GP surrogate was implemented. A GP is a random process that is characterized by sampling functions f(ξ) such that at a given point [small xi, Greek, circumflex] the samples f([small xi, Greek, circumflex]) have a normal distribution,where μ(ξ) is a mean function, here taken initially to be zero because we have no prior information, and K(ξ,ξ′) is a covariance function. The GP has 5 inputs, d, d_std, S, S_std which are the parameters that define the geometry of the microstructure, and λ which is the stretch of interest. These inputs are denoted with ξ to distinguish them from the coordinates x. The covariance matrix is based on a radial basis function (RBF) kernel following,^30,51 under the assumption of smoothness of the stress as a function of the inputs. The components of the covariance matrix between two points ξ⁽ⁱ⁾ and ξ^(j) arewith Λ = diag(Λ₁, Λ₂, Λ₃, Λ₄, Λ₅), and each Λ_i capturing the squared characteristic length-scale of each input and ζ denoting other hyperparameters of the kernel such as the signal strength.

Note that due to the assumption of the noisy observations we have an extended covariance Σ = K + σ_ms²I, which is the sum of the GP covariance matrix plus Gaussian noise σ_ms². Given a set of input data Ξ = [ξ⁽¹⁾, ξ⁽²⁾, … ξ^(N)], and corresponding observations of the stress at those inputs Y = [σ^e(1), σ^e(2), … σ^e(N)], the predictive mean and variance for a new input ξ* are

μ_pred (ξ⁽*⁾) = k^TΣ⁻¹Y, (14)

σ_pred²(ξ⁽*⁾) = k(ξ⁽*⁾, ξ⁽*⁾;ζ) +σ_ms² − k^TΣ⁻¹k. (15)

In (14), the vector k is

k = (k(ξ⁽*⁾, ξ⁽¹⁾;ζ), …, k(ξ⁽*⁾, ξ^(N); ζ)). (16)

2.5 Finite element model

We aim to replace the microstructural model of the subcutaneous tissue with an equivalent homogenized response to enable finite element simulations. This will enable us to find effective elastic properties to perform complex simulations related to drug injection physiology while using the inferred behavior due to the microstructure. To accomplish this, we assume that the subcutaneous tissue behaves as a Blatz–Ko material with strain energy defined as³²in terms of the invariants of deformation I_i, and the parameters are the shear modulus μ, the Poisson ratio ν, and the scaling parameter ϕ. The elastic stress followswhere , I is the identity matrix, b = FF^T is the left Cauchy Green deformation tensor. We can then apply deformation gradient (9) to eqn (17) and (18) to obtain the solution for the equi-triaxial problem.

We can find the optimal μ, ϕ, and v by minimizing the root mean square error

where σ^e(i) are the stresses predicted by the GP (the mean prediction μ_pred) and [small sigma, Greek, circumflex]^e(i) are the predicted stresses with the Blatz–Ko model for a given triaxial input λ and parameters μ, ϕ, and v.

The Blatz–Ko model fitted to the GP has the advantage that it generalizes stress predictions for arbitrary b and not just tri-axial deformation. It can be implemented as a user material into finite element packages such as the popular FEBio package.³¹ To model drug injections, the model has to be poroelastic and, unlike the equilibrium equations solved for the microstructure, has to account for fluid flux. Following extensive work in the literature by us and others,⁵² permeability can be effectively modeled as a function of the volume change J = det F using the Holmes–Mow model.²⁰ We assume the mean permeability parameters reported for porcine subcutaneous tissue.⁵² We perform a simulation where a subcutaneous tissue sample sees a pressure similar to an injection pressure P = 2 kPa. The material parameters used for the simulations are the stiffer and softer materials reported in Table 2. Since permeability is out of the scope of this study it will be assumed to be independent of RVE microstructure.

Table 2 Blatz–Ko effective elastic parameters

	d = 59.6 ± 7.8 μm	d = 75.2 ± 24.28 μm	d = 105.6 ± 14.7 μm
Parameter	S = 0.82 ± 0.0306	S = 0.82 ± 0.0306	S = 0.82 ± 0.0306
μ [kPa]	2.329	2.222	1.777
ϕ	0	0	0
ν	0.397	0.388	0.386
RMSE [kPa]	0.0455	0.0136	0.0253

---

3 Results

3.1 RVE deformation and homogenized response depend on packing geometry

The RVEs captured the entire range of geometries according to Table 1 and illustrated in Fig. 1. When subjected to tri-axial deformation, individual cells showed highly heterogeneous deformation. For instance, as shown in Fig. 2, for a tri-axial deformation given by λ = 1.15, the corresponding average volume change that is expected in the RVE is J = 1.52. However, as seen in Fig. 2b, for a random RVE in the set, the individual change in cell volume ranges from almost no volume change J^(c) = 1, to extreme deformations of J^(c) = 4.13. The average is indeed 〈J^(c)〉 = 1.52 as expected from the boundary conditions. The area, similarly, has a wide distribution, from J^(c)_a = 1 to J^(c)_a = 2.56 for this particular RVE (Fig. 2d). Variation in volume and area changes stems from the variability in the cell shape and size within an RVE which can lead to cell jamming.^53,54 In the extreme case of a perfect packing, all cells would deform equally. This underscores the importance of sampling random packings as occurs in adipose tissue.

Having generated 1200 RVEs spanning the geometries in Table 1 and done simulations on 1068 of them (based on the cutoff density), we can ask how exactly do the adipose tissue mechanics depend on geometry? The GP was trained on the simulation data from the RVEs at the five different stretches λ = [1.0, 1.05, 1.10, 1.15, 1.20] for a total of 2359 data points. The GP surrogate then allowed us to interpolate the stress as a function of geometry and deformation. As depicted in Fig. 3a, stress increases nonlinearly with stretch for the mean values of the other parameters. Fig. 3b and d show the dependence on cell size and sphericity (with all other inputs set at their mean value). Size and sphericity contribute significantly to the homogenized response. A decrease in cell size leads to significantly higher stresses somewhat proportional to d^−1/2. The proportionality of mechanical properties to d^−1/2 is called the Hall–Petch relationship and it was established originally for the yield stress of metals, but extended to other mechanical properties since then.^55–57 In the adipocyte system, the increase in stress with d^−1/2 is likely due to the increase in collagen area with respect to the adipocyte volume. An increase in sphericity reduces the stresses. Fig. 3g shows the interplay between sphericity and cell size and shows that cell size has a greater effect. In particular, for high cell sizes, the sphericity effect becomes almost negligible. For the lower range of cell size, on the other hand, sphericity can become important in tissue softening.

Fig. 3 Stress response by fixing all parameters but (a) λ, (b) d, (d) S. The shaded blue region in (a), (b) and (d) shows the 95% confidence interval from the GP. (c), (e)–(i) Contour plots for stress when two of the inputs are varied while keeping the rest at their mean value and λ = 1.15.

The effect of the other parameters and their pairwise couplings are shown in Fig. 3c, e–i. We show that the standard deviations of the geometry can have an effect, although less pronounced, on the homogenized response. For example, the standard deviation of the cell size d_std has a comparable importance to sphericity (Fig. 3e). Distributions of cell size with higher variance (higher d_std), increase the stress with respect to sphericity, see Fig. 3e. For metals, wider distributions of grain size (higher d_std), reduce the yield strength when compared to mean grain size d.^58,59 This is not fully evident in our system. In Fig. 3h, the size d has the dominant effect. For smaller adipocytes, increasing the width of the distribution indeed increases the stress slightly with respect to d as illustrated in Fig. 3h.

The other parameter characterizing the variance in adipocyte packing geometry is the standard deviation of the sphericity S_std. S_std has an influence on the elastic response comparable to sphericity S (Fig. 3f). For lower sphericities S, the stresses are higher. For packings with more spherical cells on average, increasing the variance of the sphericity distribution reduces the stress even further (Fig. 3f). Experimental characterization has reveled correlations between sphericity and standard deviation of sphericity for several epithelial tissues, suggesting that not all combination of geometric parameters are equally likely in nature.⁵⁴ Recent experimental work measuring stiffness of epithelial packings with atomic force microscopy (AFM) has shown that reducing sphericity indeed results in stiffer tissue.⁶⁰

The bottom line, however, is that cell size is the most important contributor to the homogenized response (Fig. 3g–i), with the other three geometric parameters having smaller but comparable sensitivity.

We generated a very wide range of microstructures and interpolated them with the GP surrogate in order to predict the stress–stretch response given a new, potentially unseen, microstructure. In particular, we can approximate geometry parameters from histology images and use our GP surrogate to estimate their equi-triaxial deformation. The reported values for statistics of different histology images for subcutaneous tissue in previous studies^8,18 provide the necessary input to predict the stress–stretch response for those histology images, as shown in Fig. 4.

Fig. 4 Equivalent diameter and width of the size distribution have an effect on the stress response. Geometric values were taken from thorough characterization of adipocyte size in the literature. A sphericity of 0.82 ± 0.0306 was assumed.¹³ The shaded regions depict the 95% confidence interval. Histology slide was taken from the experimental work⁶¹ on porcine adipose tissue to showcase the typical packing of adipocytes in this tissue.

3.2 Finite element simulations of drug delivery reflect microstructure effects

Even though we only explore the tri-axial deformation in the RVE cases, as mentioned in the Methods section, we can extend to arbitrary deformations by fitting a constitutive model such as Blatz–Ko.³² Extensions of the GP to arbitrary deformations would be possible by building GP energies in terms of invariants of deformation rather than only tri-axial stretch, and imposing suitable regularizations such as polyconvexity and other regularizations such as minimizaiton of second derivatives (to promote linear behavior outside the training region).^2,62,63 We found that the Blatz–Ko model was actually a very accurate model for tri-axial deformation as see in Fig. 5a and b. we can find the optimal Blatz–Ko material parameters for a few RVEs as shown in Table 2. The model was implemented as a user material in FEBio. Fig. 5c shows the corresponding mesh used to verify the implementation by doing the quasi-static inflation with zero traction at all faces. The uniform solution behaves exactly as the GP.

Fig. 5 (a) RVE characterized by d = 59.6 ± 7.8 μm and S = 0.82 ± 0.0306, (b) corresponding GP prediction, which aligns with the FEBio implementation of the Blatz–Ko material for the case of uniform fluid pressure loading (c).

We picked two extremes corresponding to very different RVE microstructures. For a microstructure d = 59.6 μm, d_std = 7.8 μm, S = 0.82, S_std = 0.0306, captured with a Blatz–Ko model μ = 2.329 kPa, ν = 0.397, ϕ = 0, the tissue homogenized properties are stiff. We simulated an injection pressure of 2 kPa at the center of a 1 × 1 × 1 mm³ domain with symmetric boundary conditions on the sides, fixed vertical displacement at the bottom boundary, and zero fluid pressure at the top surface as seen in Fig. 6a. The pressure leads to fluid flux (arrows in Fig. 6a, right column), which expands the tissue near the injection point. Over time, the pressure iso-contours extend outwards as the whole domain expands toward equilibrium. Volume changes at the injection site are J = 1.15, decaying to 1.11 and 1.07 by 0.22 and 0.43 mm respectively.

Fig. 6 (a) Injection into a stiff adipose tissue domain based on the packing geometry d = 59.6 μm, d_std = 7.8 μm, S = 0.82, S_std = 0.0306 corresponding to a stiff microstructure. Injection is controlled by a pressure of 2 kPa at the center of the domain, which leads to local volumetric expansion which extends to the surrounding tissue as time progresses. (b) The same boundary conditions but for a packing specified by d = 105.6 μm, d_std = 14.7 μm, S = 0.82, S_std = 0.0306 result in larger deformations and larger fluxes.

For a microstructure characterized by d = 105.6 μm, d_std = 14.7 μm, S = 0.82, S_std = 0.0306, the material is soft and captured with Blatz–Ko parameters μ = 1.777 kPa, ν = 0.386, ϕ = 0 (Fig. 6b). Injection into this subcutaneous tissue material of a fluid pressure P = 2 kPa at the center of a domain leads to volume changes up to J = 1.26 near the injection point. Because volume changes increase permeability,⁶⁴ the pressure decays rapidly away from the injection point, accompanied by larger fluid fluxes. Overall, the softer tissue acquires greater fluid volume for the same pressure. Even at 0.45 mm from the injection site at t = 10 seconds, the volume change is J = 1.2, indicating that the fluid has reached almost the entirety of the domain. It is also worth noting that the Green Lagrange strain tensor at the injection location was close to triaxial, justifying the training of the model based on triaxial deformation of the RVEs. For the stiffer material, the strain at the center was E_xx = 0.045, E_yy = 0.045, E_zz = 0.054, E_xy = 1.9 × 10⁻⁴, E_xz = 6.1 × 10⁻⁵, E_yz = −6.7 × 10⁻⁵, while for the softer material the strain at the center was E_xx = 0.067, E_yy = 0.067, E_zz = 0.079, E_xy = 1.9 × 10⁻⁴, E_xz = 2.7 × 10⁻⁴, E_yz = −4.5 × 10⁻⁴. In summary, just by changing the microstructure, the response during drug injection can be remarkably different as illustrated in Fig. 6, even if all the material parameters of individual constituents are kept fixed.

4 Discussion

In this study we constructed 1200 RVEs of Laguerre packings to describe the possible geometries of subcutaneous tissue.^28,29 Adipose tissue microstructure was parameterized by mean cell diameter, sphericity, and standard deviation these two geometric quantities. The RVEs were subjected to tri-axial loading, an expected deformation mode during drug delivery, whereby adipose tissue locally expands near the injection site due to an increase in fluid pressure.² The equilibrium configurations of the cells in the packing upon deformation were obtained by minimizing the energy of the system and imposing periodic boundary conditions. After sampling a large range of possible microstructures, we were able to the interpolate between the geometry and the homogenized stress response with the aid of a GP surrogate.⁶⁵ We also extended from the RVE simulation to finite element simulations of drug delivery leveraging the Blatz–Ko constitutive model.³²

A hallmark of metabolic disease is the increase in adipocyte size with respect to healthy individuals.⁴³ Increase in adipose tissue volume due to hypertrophy rather and hyperplasia (increase adipocyte number via cell division) does correlate with adipose tissue dysfunction.^23,25 Here we showed that the increase in adipocyte size should lead to softening of the tissue all other factors remaining constant. The proportionality of d^−1/2 is known as the Hall–Petch relationship for linear elastoplastic polycrystaline materials.⁵⁷ However, we show that size alone is not the best predictor of adipose tissue biomechanics due to the nonlinear coupling to other geometry parameters. For instance, comparison of size and sphericity shows that sphericity can become important, with more spherical cells reducing the stress compared to elongated cells.⁶⁰ In metabolic disease, the reduced proliferation of cells and their enlargement can lead to more disorganized geometry, overall reducing sphericity,¹³ which would lead to an increase in stiffness, counteracting to a small extent the effect of the larger cell size. Moreover, mean size and mean sphericity are, by themselves, not fully predictive of the biomechanics of the tissue: spread of the distribution of size and sphericity is also crucial (see Fig. 3). Even though we show the importance of these four geometric parameters, there is a dearth of detailed data of the 3D packing of adipocytes and changes associated with age and disease.^8,66

One important assumption employed in this work was that material parameters of individual constituents remained constant. This allowed us to isolate the effect of geometry alone on the homogenized response. However, it is well-known that metabolic disease affects not just adipocyte size but also triggers a pro-fibrotic response that can alter the composition, thickness, and mechanical properties of the collagen foam surrounding adipocytes.^67,68 Adipocytes themselves can exhibit changes in mechanical properties as a function of lipid droplet numbers per adipocyte, a sign of adipocyte differentiation.¹⁰ Thus, future work should further extend the present model to sample the space of material properties on the homogenized response.

Applications of interest, such as drug delivery to the subcutaneous tissue, require large scale simulations for which modeling individual adipocytes would not be feasible.⁶⁹ Thus, there is a need to obtain homogenized subcutaneous tissue properties.^4,5,15,70,71 Multiscale models to tie the cell to the tissue scale are often too expensive for many realistic simulations.⁷² While it is possible to embed RVEs per integration point of a finite element mesh,^73–75 the advances in machine learning and data-driven frameworks have resulted in efficient homogenization tools such as Gaussian process surrogates used here.⁶⁵ To build a complete model at the tissue scale, the RVEs generated for this study would have to be subjected to a large range of deformations and not just tri-axial loading.^63,76 Our consideration of tri-axial stretching is motivated by the application of interest, namely subcutaneous injection.⁷ For the tri-axial deformation, we found that the Blatz–Ko model offered an almost perfect fit to the GP interpolation, and we used this constitutive model to extend the response to arbitrary loading. The Blatz–Ko model was implemented in the finite element software FEBio through a user material. Other work in the field of adipose tissue biomechanics has shown that strain energies such as the Ogden model might be needed to generalize to other deformation modes.⁷⁷ During drug injections to subcutaneous tissue, the elastic properties are only partially descriptive of the tissue, the transport properties being the other set of parameters needed.⁷⁸ Permeability of adipose tissue has been reported across a wide range.⁷⁹ Only very recent studies have narrowed down the possible permeability values for adipose tissue,⁵² but these studies have ignored any possible coupling with geometry at the microscale. Here we took the permeability parameters from the literature,⁵² but additional, non-equilibrium simulations of the RVEs are needed to tie the geometry of the packing to both the elastic as well as the transport response of the tissue.

This study is not without limitations, several of which have already been discussed, such as the emphasis on the elastic response. Other nonlinearities were ignored, such as other material models,⁷⁷ or dissipative mechanisms (viscoelasticity, damage).⁷¹ Finally, experimental work is needed to concurrently image the microscale during loading in order to establish fully calibrated multiscale models.^74,75 Additional experimental work is also urgently needed to characterize the sphericity distribution in adipose tissue and how this distribution changes with age, sex and disease.

5 Conclusions

We used a GP surrogate to interpolate the response of adipose tissue under triaxial loading across all possible 3D packings of adipocytes parameterized as Laguerre tessellations. We showed that mean diameter was the most important determinant of homogenized tissue response, with proportionality close to the well-known Hall–Petch relationship of d^−1/2. This alone is of high relevance for the clinical setting because metabolic disease leads to adipocyte hypertrophy. Moreover, the emergence of new weight-loss and type-2 diabetes treatments such as Mounjaro are administered through subcutaneous injections. We show that coupling adipocyte hypertrophy to multiscale finite element simulations leads to significant changes in drug transport. Thus, we anticipate that the present work will spark new studies to better characterize 3D adipose tissue structure and constituent properties, as well as contribute to the development of patient-specific drug delivery strategies based on individual's adipose tissue microstructure.

Author contributions

JBM conceptualized the study and established methodology, coded software, performed data curation, formal analysis, visualization and writing. LS contributed funding acquisition, supervision and writing. ABT conceptualized the study and contributed with funding acquisition, supervision, methodology, formal analysis, visualization, and writing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by Eli Lilly and Company through the Purdue-Lilly partnership. The authors have no conflict of interest to declare.

Notes and references

1. K. Comley and N. A. Fleck, Int. J. Solids Struct., 2010, 47, 2982–2990.
2. Y. Leng, A. M. Ardekani and H. Gomez, J. Mech. Phys. Solids, 2021, 155, 104537.
3. A. Kivitz, S. Cohen, J. E. Dowd, W. Edwards, S. Thakker, F. R. Wellborne, C. L. Renz and O. G. Segurado, Clin. Ther., 2006, 28, 1619–1629.
4. J. L. Calvo-Gallego, J. Domnguez, T. G. Ca, G. G. Ciriza and J. Martnez-Reina, J. Mech. Behav. Biomed. Mater., 2018, 80, 293–302.
5. Z. Sun, S.-H. Lee, B. D. Gepner, J. Rigby, J. J. Hallman and J. R. Kerrigan, J. Mech. Behav. Biomed. Mater., 2021, 113, 104112.
6. M. A. Biot, J. Appl. Phys., 1955, 26, 182–185.
7. M. de Lucio, Y. Leng, A. Hans, I. Bilionis, M. Brindise, A. M. Ardekani, P. P. Vlachos and H. Gomez, J. Mech. Behav. Biomed. Mater., 2023, 138, 105602.
8. J. Honecker, D. Weidlich, S. Heisz, C. M. Lindgren, D. C. Karampinos, M. Claussnitzer and H. Hauner, Int. J. Obes., 2021, 45, 2108–2117.
9. N. Shoham, P. Girshovitz, R. Katzengold, N. T. Shaked, D. Benayahu and A. Gefen, Biophys. J., 2014, 106, 1421–1431.
10. M. Ben-Or Frank, N. Shoham, D. Benayahu and A. Gefen, Biomech. Model. Mechanobiol., 2015, 14, 15–28.
11. T. L. Alderete, F. R. Sattler, X. Sheng, J. Tucci, S. D. Mittelman, E. G. Grant and M. I. Goran, Int. J. Obes., 2015, 39, 183–186.
12. A. Michaud, J. Tordjman, M. Pelletier, Y. Liu, S. Laforest, S. Noël, G. Le Naour, C. Bouchard, K. Clément and A. Tchernof, Int. J. Obes., 2016, 40, 1823–1831.
13. J. Geng, X. Zhang, S. Prabhu, S. H. Shahoei, E. R. Nelson, K. S. Swanson, M. A. Anastasio and A. M. Smith, Sci. Adv., 2021, 7, eabe2480.
14. Y. Leng, M. de Lucio and H. Gomez, Comput. Methods Appl. Mech. Eng., 2021, 384, 113919.
15. G. Sommer, M. Eder, L. Kovacs, H. Pathak, L. Bonitz, C. Mueller, P. Regitnig and G. A. Holzapfel, Acta Biomater., 2013, 9, 9036–9048.
16. V. D. Sree, J. D. Toaquiza-Tubon, J. Payne, L. Solorio and A. B. Tepole, Ann. Biomed. Eng., 2023, 1–14.
17. C. A. Ibáñez, M. Vázquez-Martnez, J. C. León-Contreras, L. A. Reyes-Castro, G. L. Rodrguez-González, C. J. Bautista, P. W. Nathanielsz and E. Zambrano, Front. Physiol., 2018, 9, 1571.
18. V. A. Palomäki, V. Koivukangas, S. Meriläinen, P. Lehenkari and T. J. Karttunen, Adipocyte, 2022, 11, 99–107.
19. J. M. Carcione, C. Morency and J. E. Santos, Geophysics, 2010, 75, 75A229–75A243.
20. G. A. Ateshian, W. Warden, J. Kim, R. Grelsamer and V. C. Mow, J. Biomech., 1997, 30, 1157–1164.
21. K. Comley and N. Fleck, Int. J. Impact Eng., 2012, 46, 1–10.
22. S. Laforest, J. Labrecque, A. Michaud, K. Cianflone and A. Tchernof, Crit. Rev. Clin. Lab. Sci., 2015, 52, 301–313.
23. D. Lemonnier, et al. , J. Clin. Invest., 1972, 51, 2907–2915.
24. C.-Y. Lin, G. P. Sugerman, S. Kakaletsis, W. D. Meador, A. T. Buganza and M. K. Rausch, bioRxiv, 2023, 2023-03.
25. M. Varghese, J. Song and K. Singer, Mech. Ageing Dev., 2021, 199, 111563.
26. E. Chang, M. Varghese and K. Singer, Curr. Diabetes Rep., 2018, 18, 1–10.
27. J. Esque, M. S. Sansom, M. Baaden and C. Oguey, Sci. Rep., 2018, 8, 13540.
28. A. Wouterse, S. R. Williams and A. P. Philipse, J. Phys.: Condens. Matter, 2007, 19, 406215.
29. A. V. Kyrylyuk and A. P. Philipse, Phys. Status Solidi A, 2011, 208, 2299–2302.
30. C. Stowers, T. Lee, I. Bilionis, A. K. Gosain and A. B. Tepole, J. Mech. Behav. Biomed. Mater., 2021, 118, 104340.
31. S. A. Maas, B. J. Ellis, G. A. Ateshian and J. A. Weiss, J. Biomech. Eng., 2012, 134, 1.
32. P. J. Blatz and W. L. Ko, Trans. Soc. Rheol., 1962, 6, 223–252.
33. V. Sree, X. Zhong, I. Bilionis, A. Ardekani and A. B. Tepole, J. Mech. Behav. Biomed. Mater., 2023, 140, 105695.
34. P. Gómez-Gálvez, P. Vicente-Munuera, S. Anbari, J. Buceta and L. M. Escudero, Development, 2021, 148, dev195669.
35. S. Grosser, J. Lippoldt, L. Oswald, M. Merkel, D. M. Sussman, F. Renner, P. Gottheil, E. W. Morawetz, T. Fuhs and X. Xie, et al. , Phys. Rev. X, 2021, 11, 011033.
36. S. Kaliman, C. Jayachandran, F. Rehfeldt and A.-S. Smith, Front. Physiol., 2016, 7, 551.
37. A. Spettl, T. Werz, C. E. Krill and V. Schmidt, J. Stat. Phys., 2014, 154, 913–928.
38. S. T. Abraham and S. S. Bhat, Mater. Sci. Eng., A, 2023, 877, 145155.
39. T. Wejrzanowski, J. Skibinski, J. Szumbarski and K. Kurzydlowski, Comput. Mater. Sci., 2013, 67, 216–221.
40. R. Quey and L. Renversade, Comput. Methods Appl. Mech. Eng., 2018, 330, 308–333.
41. K. Lu, S. Xie, S. Han, J. Zhang, X. Chang, J. Chao, Q. Huang, Q. Yuan, H. Lin and L. Xu, et al. , J. Transl. Med., 2014, 12, 1–14.
42. C.-C. Shih, C.-H. Lin, Y.-J. Lin and J.-B. Wu, J. Evidence-Based Complementary Altern. Med., 2013, 2013, 1.
43. J. J. Fuster, N. Ouchi, N. Gokce and K. Walsh, Circ. Res., 2016, 118, 1786–1807.
44. R. Quey, P. Dawson and F. Barbe, Comput. Methods Appl. Mech. Eng., 2011, 200, 1729–1745.
45. N. Slomka, C. W. Oomens and A. Gefen, J. Mech. Behav. Biomed. Mater., 2011, 4, 1559–1566.
46. O. C. Zienkiewicz and R. L. Taylor, The finite element method for solid and structural mechanics, Elsevier, 2005.
47. G. A. Holzapfel, Meccanica, 2002, 37, 489–490.
48. J. R. H. Otter, A. C. Cassell and R. E. Hobbs, Proc. - Inst. Civ. Eng., 1966, 35, 633–656.
49. C. B. Barber, D. P. Dobkin and H. Huhdanpaa, ACM Trans. Mathematical Software (TOMS), 1996, 22, 469–483.
50. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt and SciPy 1.0 Contributors, Nat. Methods, 2020, 17, 261–272.
51. C. M. Bishop and N. M. Nasrabadi, Pattern recognition and machine learning, Springer, 2006, vol. 4.
52. P. Shrestha and B. Stoeber, Phys. Fluids, 2020, 32, 011905.
53. L. Atia, J. J. Fredberg, N. S. Gov and A. F. Pegoraro, Cells Dev., 2021, 168, 203727.
54. L. Atia, D. Bi, Y. Sharma, J. A. Mitchel, B. Gweon, S. A. Koehler, S. J. DeCamp, B. Lan, J. H. Kim and R. Hirsch, et al. , Nat. Phys., 2018, 14, 613–620.
55. E. Hall, Proc. Phys. Soc., London, Sect. B, 1951, 64, 747.
56. N. Petch, J. Iron Steel Inst., 1953, 174, 25–28.
57. N. Hansen, Scr. Mater., 2004, 51, 801–806.
58. B. Raeisinia, C. Sinclair, W. Poole and C. Tomé, Modell. Simul. Mater. Sci. Eng., 2008, 16, 025001.
59. P. Lehto, H. Remes, T. Saukkonen, H. Hänninen and J. Romanoff, Mater. Sci. Eng., A, 2014, 592, 28–39.
60. X. Xie, F. Sauer, S. Grosser, J. Lippoldt, E. Warmt, A. Das, D. Bi, T. Fuhs and J. A. Käs, Soft Matter, 2024, 20, 1996–2007.
61. M. C. Brindise, K. Buno, L. Solorio and P. P. Vlachos, Ann. Biomed. Eng., 2023, 51, 443–455.
62. N. N. Vlassis, R. Ma and W. Sun, Comput. Methods Appl. Mech. Eng., 2020, 371, 113299.
63. J. N. Fuhg, M. Marino and N. Bouklas, Comput. Methods Appl. Mech. Eng., 2022, 388, 114217.
64. V. C. Mow, S. Kuei, W. M. Lai and C. G. Armstrong, J. Biomech. Eng., 1980, 102, 73–84.
65. Y. Leng, V. Tac, S. Calve and A. B. Tepole, Comput. Methods Appl. Mech. Eng., 2021, 387, 114160.
66. S. D. Parlee, S. I. Lentz, H. Mori and O. A. MacDougald, Methods in enzymology, Elsevier, 2014, vol. 537, pp. 93–122.
67. G. Anvari and E. Bellas, Sci. Rep., 2021, 11, 21473.
68. N. Di Caprio and E. Bellas, Adv. Biosyst., 2020, 4, 1900286.
69. K. M. Moerman, M. van Vijven, L. R. Solis, E. E. van Haaften, A. C. Loenen, V. K. Mushahwar and C. W. Oomens, Comput. Methods Biomech. Biomed. Eng., 2017, 20, 483–491.
70. N. Alkhouli, J. Mansfield, E. Green, J. Bell, B. Knight, N. Liversedge, J. C. Tham, R. Welbourn, A. C. Shore and K. Kos, et al. , Am. J. Physiol.: Endocrinol. Metab., 2013, 305, E1427–E1435.
71. A. Hatt, R. Lloyd, B. Bolsterlee and L. E. Bilston, J. Mech. Behav. Biomed. Mater., 2023, 105924.
72. N. Shoham, A. Levy, N. Shabshin, D. Benayahu and A. Gefen, Biomech. Model. Mechanobiol., 2017, 16, 275–295.
73. T. Stylianopoulos and V. H. Barocas, 2007.
74. E. A. Sander, T. Stylianopoulos, R. T. Tranquillo and V. H. Barocas, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 17675–17680.
75. N. J. Witt, A. E. Woessner, K. P. Quinn and E. A. Sander, J. Biomech. Eng., 2022, 144, 041008.
76. Y. Heider, K. Wang and W. Sun, Comput. Methods Appl. Mech. Eng., 2020, 363, 112875.
77. L. A. Mihai, L. Chin, P. A. Janmey and A. Goriely, J. R.Soc., Interface, 2015, 12, 20150486.
78. A. Wahlsten, M. Pensalfini, A. Stracuzzi, G. Restivo, R. Hopf and E. Mazza, Biomech. Model. Mechanobiol., 2019, 18, 1079–1093.
79. Z. Sun, B. D. Gepner, P. S. Cottler, S.-H. Lee and J. R. Kerrigan, J. Biomech. Eng., 2021, 143, 070803.

---