Comput Mech DOI 10.1007/s00466-010-0551-8
ORIGINAL PAPER
A novel technique to parametrize shell-like deformations inside biological membranes
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
25 26 27
Abstract Biological structures like organs or embryos usually have irregular forms. The analytical description of their geometry constitutes then a major issue when numerical modeling is used to simulate their behavior. The novel technique presented here allows to parametrize many threedimensional objects whose shape can be approximated by the one of a general shell. The main idea of the present work lies on the coupling between electricity, geometry and mechanics because the electric potential inside an object has properties inherited from the geometrical representation. Based on this fundamental property, we are able to describe both simply and multiply connected domains. Two forms of the theoretical formulation are proposed. First, we present a partial version for which the electric potential is employed to compute only one of the principal variables, which compose the general curvilinear system. Second, we introduce an extended version for which the complete covariant and contravariant bases are obtained and the conventional kinematics and mechanics for a thick general shell are redefined. Finally, for illustrative purpose, we propose a biomechanical application on a rather complex structure such as Drosophila embryo. In particular, we show how the new technique fits well when a parametrical description of the system geometry is required to simulate the morphogenetic movements.
1 Introduction
28
Very often in life sciences we find complex and irregular forms. The analytical description of their geometry represents a major issue when computer modeling is used to simulate their behavior. Furthermore, many biological systems such as tissues and membranes may deform themselves through a series of successive processes, so that their shape becomes more and more complicated. Actually, only Mother Nature knows the initial computeraided design (CAD) representation of any living organism, whereas for computer simulations good approximations of the real geometrical shapes are usually achievable. Besides, in the last decades, recent new imaging techniques such as magnetic resonance imaging (MRI) and constant time imaging (CTI) together with several microscopy methods have played a significant role to provide digital images of living tissues, even though many parameters are still overlooked or unknown. There are some softwares which can supply a sufficiently accurate finite element mesh for several anatomical (i.e. the heart, the scapula, the ankle, the knee or the brain) or biological structures (i.e. the cell or the embryo). The main drawback of this approach is that a parametrical study on the influence of the shape and other geometrical properties such as the thickness or the curvature is often very difficult to carry out. Although they vary from one case to the other, having these parameters available, even with an approximate geometry, would certainly allow us better understanding the main role of each sub-structure. Furthermore, very often several structural components may have either one small dimension (i.e. membranes like shells) or two (i.e. fibers like networks). Thus, in order to have a precise reproduction of the strains such as bending or twisting, the middle surface in the first case or the middle line in the second case need to be modeled. Shell-like or
cted
2
orre
1
unc
Author Proof
Received: 5 July 2010 / Accepted: 25 October 2010 © Springer-Verlag 2010
pro of
R. Allena · D. Aubry
Keywords Biological shells · Electric potential · Simply and multiply connected domains · Morphogenesis
R. Allena (B) · D. Aubry Laboratoire MSSMat, Ecole Centrale Paris, Grande Voie des Vignes, 92295 Châtenay-Malabry, France e-mail:
[email protected]
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
TYPESET
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
Comput Mech
66 67 68 69 70 71 72 73 74 75 76 77
Author Proof
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
138
2 Review of basic principles for shell mechanics
139
pro of
65
This is actually the key idea of the present work: to parametrize one of the coordinates of a general shell, specifically the thickness variable, using the electric potential and to define then a curvilinear coordinates system, which enables to describe any type of three-dimensional geometry, even those showing irregularities or holes. Eventually, an extension of the technique will even allow building tangential curvilinear coordinates. Besides, the approach provides an analytical formulation of kinematic quantities very much reduced and more intuitive compared to the one adopted in our previous work [3]. Finally, we will show here that the applications and the relative results are still very relevant and conclusive with respect to previous works and the reality. The paper is composed of four principal sections. In the first one, we recall the general concepts needed for shell theory through an example of a regular geometry. The second section is focused on the detailed description of the new approach. Two forms of the formulation are proposed: a reduced and a complete version. In the forth part of the work, the new capacitive electrical technique is used to propose an alternative geometrical parametrization and to describe the mechanics of a general shell. Finally, in order to illustrate the method, we present a biomechanical application on a rather complex structure such as the Drosophila embryo (Table 1).
cted
64
rod-like descriptions are then mandatory and digital segmentation usually does not directly provide such a representation. Indeed, modeling those membranes as shells is very interesting apart from the reduction of the number of degrees of freedom. Actually, biologists have often identified, at least qualitatively, elementary deformations typical of those thin structures. For instance, they have been able to observe convergence–extension movements, apical constrictions or apico-basal elongations of embryonic tissues during morphogenesis [15,18], which can be translated into a superposition of shell elementary deformations mentioned above. Obviously, this modeling may be achieved on highly idealized structures such as spheres [2], ellipsoids [3] or other simple analytical geometries. But real biological geometries have frequently non-uniform thickness and curvatures, which implies that this mathematical approach, although very interesting for the understanding of basic principles, is in our opinion too limited. Consequently, we would like to develop an intermediary technique, which lies at the middle line between an analytical approach for simplified biological structures and a digital segmentation of real complex living structures. From a structural mechanics standpoint, a shell is a geometrical body endowed with several features such as a middle surface and two outer surfaces so that a normal to the former intersects the outer surfaces at a distance h/2, not necessarily constant [20]. Then, a dimensional reduction of the number of d.o.f is usually obtained by substituting appropriate averages of the dependent variables with respect to the fiber variable. Building reasonable middle surfaces or lines of complex membranes or fibers, respectively, and computing the elementary deformations acting on them is thus considered as an important tool to increase the link between mechanics and biology. Here, we propose an original approach based on the properties of the electrostatic potential, which is known having a strong connection with the geometry of three-dimensional objects [16], even the most complex, inside which it is calculated. This method can be very useful for two main reasons. First, the coupling of mechanics and electrics may be compulsory when trying to model structures such as nerves or muscles for which the electrical activity may play a significant role. Thus, be prepared for it would allow taking into account potential electrical signals also at the interior of other biological systems. So far this aspect has been studied and numerically simulated at least for the case of the cardiac tissues [13]. For the embryo, and more specifically for the cell, this issue is still little known and explored. Second, as already mentioned, the electric potential has been often recognized as a quantitative witness of the geometrical shape. For instance, every college student certainly remembers the relationship in physics between the solid angle and the potential of a unit charge through a closed surface or between the volume of a conductor and its electric capacity.
Let us consider a thick three-dimensional membrane 0 as represented in Fig. 1a. In this particular case, one dimension—the thickness—is about an order of magnitude smaller than the others (i.e. the diameter). The membrane may be
orre
63
unc
62
Table 1 Main parameters of the model Variable
Notation Default value
Young’s modulus
E
100 Pa
Poisson’s ratio
ν
0.45
Lamé’s coefficient
λ
310.34 Pa
Lamé’s coefficient
µ
34.48 Pa
Ellipsoid semi axis
a
250 µm
Ellipsoid semi axis
b
75 µm
Ellipsoid semi axis
c
75 µm
Ellipsoid thickness
h
15 µm
Refined ellipsoid major axis –
500 µm
Refined ellipsoid semi axis CE
175 µm
Refined ellipsoid semi axis DF
165 µm
Active region VFI
–
Amplitude along Vθ 50◦
Active region GBE
–
Length along AP axis 350 µm Length along AP axis 200 µm Amplitude Vθ 90◦
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
TYPESET
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
140 141 142 143
Comput Mech
Author Proof
148 149 150 151 152 153 154 155 156
157
158 159 160 161 162 163 164
165
166 167 168 169
170
171 172 173 174 175 176 177
178
p = p0 (ξ1 , ξ2 ) + ξ3 n0
form now the covariant basis at x. Thus, Eq. (5) can be now rewritten as
(1)
where p0 represents any point on the middle surface 0 of the membrane, ξ3 is the thickness of the membrane and n0 is the normal vector to the middle surface passing by p0 . This parametrization of the shell leads to a system of curvilinear ∂p , the partial coordinates with the covariant basis gp,m = ∂ξ m derivative of p with respect to ξm , and the contravariant basis gm p = ∇p ξm , the gradient of ξm with respect to p. It can be shown that, using the Einstein summation convention, we have
gx,m =
179
∂x ∂ξm
(6)
F = g x,m ⊗ g mp
(7)
With the covariant basis at x and the contravariant basis at p, we are able to describe any type of deformation. The deformation gradient components Fn,m can be obtained taking the projection of the deformed vectors gx,m onto the undeformed vectors gp,n , so that Eq. (7) becomes n m (8) F = gx,m , gp,n gnp ⊗ gm p = Fn,m gp ⊗ gp
cted
147
where
gp,m ⊗ gm p =I
(2)
where I is the identity tensor and the tensor product between two vectors a and b is defined as follows (a ⊗ b)i j = ai b j . Without inverting the initial configuration in the real space domain p from the one in the parameter domain ξ , we would like to find the vectors gm p . According to this idea, the contravariant basis is not deduced from Eq. (2), but from the orthogonality property, which reads gm (3) p , gp,n = δmn
where gnp ⊗ gm p is the initial basis as defined above. This approach fits well as long as the transformation is simple and Fn,m can be given a priori. To conclude, let us consider the particular case of a hollow ellipsoid. This geometry will also be employed later on in this paper for an interesting biomechanical application to the Drosophila embryo (Sect. 5), as similarly as in our previous work [3], and will eventually illustrate our new technique. Using a cylindrico-polar basis (r, θ, z) where
orre
146
open or closed, simply or multiply connected. Using a curvilinear coordinates system (ξ1 , ξ2 , ξ3 ), any point p at its reference configuration is defined [20] as
where the scalar product for two vectors a and b is defined as follows (a, b) = i ai bi . Therefore, through a series of mathematical manipulations [11], we obtain i.e. gp,2 ∧ gp,3 (4) g1p = gp with gp = gp,1 ∧ gp,2 , gp,3 . The other two vectors g2p and g3p composing the contravariant basis are obtained by simple index permutation. When the shell undergoes a transformation from 0 to , any point p deforms so that it comes to its present position x and the relative deformation tensor F is classically expressed as
unc
144 145
pro of
Fig. 1 a Parametrization of a general shell. b Geometrical representation of a regular ellipsoid using a curvilinear coordinates system
∂x F= ⊗ ∇ p ξm ∂ξm
(5)
ir = (cos θ, sin θ, 0)
iθ = (− sin θ, cos θ, 0)
(9)
iz = (0, 0, 1)
p0 (θ, z) = ρ (z) ir (θ ) + ziz
(10)
with a and b the semi axes of the ellipsoid. Since some deformations may be described with respect to the middle surface of the structure, it is necessary to introduce the reference position of any point p through the thickness of
LE
183
184 185 186 187 188
189
190 191 192 193 194 195 196 197 198
200
202
where ρ (z) parametrically describes the geometry and is equal to z 2 (11) ρ (z) = b 1 − a
DISK
182
201
any point p0 on the middle surface 0 is expressed as
TYPESET
181
199
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
180
CP Disp.:2010/11/4 Pages: 16 Layout: Large
203
204 205
206
207 208 209 210
Comput Mech
212
the ellipsoid as p (θ, z, ζ ) = p0 (θ, z) + ζ n0 (θ, z) = ρ (z) ir (θ ) + ziz + ζ
214
215 216 217 218
219
220
221
−ρ ′ (z) iz + ir (θ ) √ 1 + ρ ′′ (z)
with
224
225 226 227
228
229
230
231 232 233 234 235 236 237
aθ = ρ +
1 + (ρ ′ )2
ζρ ′′ az = 1 − 3 1 + (ρ ′ )2 2 aζ = 1 + ρ ′′
with α ben once more controlling the amplitude of the applied strain. In Fig. 2b, the results show the twisting of the structure depending on the thickness of the ellipsoid. Here, we have reviewed the basic principles to describe the geometry and the kinematics of thin shells knowing a priori the relative curvilinear coordinates system. In the next section, we are going to propose a new method, which allows building a special coordinates system relevant even for those more complex forms and membranes.
238 239 240 241
242
243 244
245
246 247 248 249 250 251 252 253 254
(14)
3 Capacitive electric geometry
Let us imagine now applying to the ellipsoid a uniform but curvilinear elongation in the z direction. The final position x can then be written as ′ −ρ (z) iz + ir (θ ) (15) x (θ, z, ζ ) = p0 (θ, z) + ζ √ 1 + ρ ′′ (z) where z = 1 + α el z
where the new angle θ is now function of both θ and ζ as follows θ = 1 + α ben ζ θ (18)
orre
223
(12)
where ζ is the distance between p and its projection p0 on the middle surface (Fig. 1b). Then, the deformation gradient F will be expressed, according to Eq. (5), as ′ ρ ir + iz 1 ∂x 1 ∂x ⊗ iθ + ⊗ F= aθ ∂θ az ∂z 1 + (ρ ′ )2 1 ∂x ′ ⊗ −ρ i z + i r (13) + aζ ∂ζ ζ
222
unc
Author Proof
213
also because describing it in a curvilinear system of coordinates is not so easy. However, our analytical formulation allows computing it rather simply. Thus, the final position x becomes equal this time to ′ −ρ (z) iz + ir θ (17) x (θ, z, ζ ) = p0 θ, z + ζ √ 1 + ρ ′′ (z)
cted
211
pro of
Fig. 2 Examples of elementary deformations on a regular ellipsoid. a Cross-section of the ellipsoid along the x-axis when an elongation in the same direction is applied. b Cross-section along the y direction when a bending is introduced
(16)
with α el controlling the amplitude of the deformation. Then, the deformation gradient F is computed according to Eq. (13). In Fig. 2a, we can observe the results for a quarter of the ellipsoid that show actually a lengthening along the horizontal axis of symmetry. The bending along the small crosssection of the ellipsoid, with the bending axis parallel to the z axis, constitutes another interesting example of deformation,
It is well known that there is a intricate connection between solutions of Laplace’s equation on a given domain with appropriate boundary conditions and the proper geometry of the same domain. Simple examples are given by the definition of the solid angle, the geometric capacity measure or the conformal mapping in two dimensions [16]. Based on these remarks, we propose here an alternate method to build the previous curvilinear bases by exploiting this basic property. This will allow us describing typical shell deformations and mechanics for those biological structures more complex than the ones that are classically represented by analytical expressions.
267
3.1 The normal coordinate electric potential
268
We restrict ourselves to a thick membrane formed by an outer boundary ∂e and an inner boundary ∂i , these two surfaces
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
255
TYPESET
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
256 257 258 259 260 261 262 263 264 265 266
269 270
Comput Mech
pro of
Fig. 3 On the left, numerical computation of the electric potential Vξ3 through the thickness of a general shell. From blue (−h/2) to red (h/2) the strength of the computed electric potential. On the right, representation of the external (red), the middle (green) and the internal (blue) surface of the shell
273 274 275 276 277 278 279 280 281 282
283
284 285 286 287 288 289 290 291
cted
272
being rather close to each other with respect to a typical overall length scale of the membrane and a lateral boundary. Let h be an average estimation of the thickness, we would like to determine the middle surface 0 and the thickness ξ3 of the membrane using as a distance parameter an electric potential with uniform boundary conditions. If div is the divergence and ∇ is the gradient, then we define the electric potential Vξ3 (x, y, z) (Fig. 3) within the thickness of the membrane through the Laplace’s equation and the Dirichlet’s boundary conditions of a dielectric material connected with electrodes at two opposite imposed potentials ⎧ −div c∇Vξ3 = 0 inside 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ Vξ3 = + h2 on the external boundary ∂e (19) h ⎪ ⎪ on the internal boundary ∂ V = − i ξ ⎪ 3 2 ⎪ ⎪ ⎪ ⎪ ⎩ ∂ Vξ3 ∂n = 0 on the lateral boundaries
approximation of the middle surface of the membrane as it can be observed in Fig. 3. Then, from the calculated electric potential Vξ3 , an approximate normal vector n 0 to the shell middle surface Vξ3 = 0 (Fig. 4a) can be computed as follows ∇Vξ3 n0 ∼ = ∇Vξ 3
(20)
orre
271
unc
Author Proof
Fig. 4 Numerical computation of the curvilinear basis. a The normal vector n0 . b, c The tangential vectors Vθ and Vz , respectively
where n is the geometrical normal vector and the coefficient c, which may even be a matrix, is assumed to be here always equal to 1. Because of the classical properties of the Laplace’s equation, it is easily checked that the isosurfaces Vξ3 = const do not cross the inner and the outer boundaries due to the maximum principle [10] and they are also parallel surfaces. Furthermore, the isosurface Vξ3 = 0 usually gives a good
Using an asymptotic expansion with respect to the thickness, it can be shown that the thinner the membrane, the more ∇Vξ3 is perpendicular to the middle surface of the shell. Consequently, Vξ3 (x, y, z) appears to be a valid candidate to replace the usual normal curvilinear coordinate ξ3 (x, y, z). Even though we do not explicitly seek here the two other tangent variables, we have to define the tangent vectors to the equipotential surfaces. In the same spirit as in the finite elements approach for thick shells formerly proposed by Zienkiewicz [27], by a simple set of cross products, we find a new curvilinear basis constituted by the tangent vectors V 1 and V2 (Fig. 4b, c) and the normal n0 . V 1 and V 2 are expressed, from the global unit vector iz and i y (or any other vector not parallel to ∇Vξ3 ) as follows i z ∧ ∇Vξ3 V1 = i z ∧ ∇Vξ 3 i y ∧ ∇Vξ3 V2 = i y ∧ ∇Vξ
TYPESET
LE
295 296
297
298 299 300 301 302 303 304 305 306 307 308 309 310 311
312
(22)
313
3
DISK
294
(21)
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
292 293
CP Disp.:2010/11/4 Pages: 16 Layout: Large
Comput Mech
315 316 317
319 320 321 322 323 324 325 326 327 328 329 330 331
g1p =
V 2 ∧ n0 (V 2 ∧ n0 , V 1 )
(23)
The associated curvilinear coordinates ξm could be recovered, at least theoretically, by integrating the differential system given by the definition of the basis. As shown in Sect. 2, when elementary transformations are introduced, the deformation gradient F can be calculated from Eq. (8). Let us apply the technique to the regular ellipsoid presented in Sect. 2. The Laplace’s equation, which supplies the electric potential Vξ3 , is numerically solved by a standard finite element method, which will provide an approximation to this unknown and to the basis. If we consider the two basic transformations presented in Sect. 2, the total deformation gradient F can be directly deduced, so that we have respectively
332
F = I + α el V 2 ⊗ V 2
333
for the elongation in the direction of V 2 and
335 336 337 338 339 340 341 342 343 344 345 346 347 348
F = I + α ben Vξ3 V 1 ⊗ V 1
3.2 A complete curvilinear coordinates system for simply and multiply connected domains When a full displacement transformation is required, an explicit system of coordinates is often necessary. We need then to build two other coordinates, complementary to Vξ3 . First, for a simply connected object (without holes), the method presented above can be easily extended to outer faces, so that the two completing variables Vξ1 and Vξ2 are numerically computed as similarly as for the electric potential through the thickness of the structure [Eq. (19)]. The corresponding Laplace’s equation and the normalized Dirichlet’s boundary conditions on the suitable lateral faces are then given by ⎧ −div c∇Vξ1 = 0 inside 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ Vξ1 = +1 on ∂1+
351 352 353 354 355 356 357 358 359 360 361
362
⎧ −div c∇Vξ2 = 0 inside 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ Vξ = +1 on ∂2+ 2
(25)
for the bending along the cross section of the shell. The numerical results are very close to the ones obtained for the ellipsoid defined by analytical coordinates (Fig. 2a, b, respectively). In Fig. 5, we have envisaged more complex membranes with curvature variations even of the order of the thickness itself, which usually constitutes a challenge for shell mechanics but are not rare in biology. Only in this extreme case, we see that the isosurface V ξ3 = 0 is somewhat drifted with respect to the middle surface. Therefore, when only one of the curvilinear coordinates is explicitly computed, a formulation is obtained, which works well for both elementary and more complicated examples, even for closed membranes.
⎪ Vξ2 = −1 on ∂2− ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ Vξ2 ∂n = 0 on the other boundaries
(26)
In Fig. 6 an example with the isovalues of the three variables is represented. Secondly, in a domain with holes, like a hollow ellipsoid (Sects. 2, 5) or even an annular domain, there are often only two to four boundaries, but not six as necessary. Thus, all the tangential variables cannot be simply calculated [6,12]. To circumvent this issue, we propose an original technique to recover the same properties for a multiply connected domain like a closed thick membrane. We restrict ourselves again to the particular hollow ellipsoid already considered. Only one of the two tangential variables is calculated, specifically θ , which will be indicated as
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
349 350
⎪ Vξ1 = −1 on ∂1− ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ Vξ1 ∂n = 0 on the other boundaries
orre
334
(24)
unc
Author Proof
318
Therefore, using Eq. (3), if we assume that gp,1 = V 1 , gp,2 = V 2 and gp,3 = n0 , which compose then the new covariant basis, we are able to find also the new contravariant basis gmp and we obtain
cted
314
pro of
Fig. 5 Influence of the curvature thickness ratio on the isovalue surface Vξ3 = 0
TYPESET
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
363
364 365 366 367 368 369 370 371 372 373 374 375
Comput Mech
pro of
Fig. 6 Isovalues of the three potentials giving the coordinates system constituted by Vξ1 (a), Vξ2 (b) and Vξ3 (c) for a simply connected membrane
Author Proof
The same technique can be further extended to calculate Vϕ , which would indicate the coordinate ϕ in the case of a spherical coordinates system. The relative Laplace’s equation in the case of a half hollow ellipsoid would be equal to
(28)
cted
⎧ −div c∇Vϕ = 0 inside 0 ⎪ ⎪ ⎨ V = 0 on ∂ ϕ ϕ0 V = π on ∂ ⎪ ϕ ϕπ ⎪ ⎩ ∂ Vϕ ∂n = 0 on the other boundaries
Fig. 7 For sake of clarity, only half of the hollow ellipsoid is reproduced. ∂s1 represents the slit surface needed for the computation of Vθ . ∂ϕ0 and ∂ϕπ are instead the two slit surfaces introduced to calculate Vϕ in the next paragraph
378 379 380 381 382 383 384 385 386 387
388
389 390 391
orre
377
Vθ . The last cylindrico-polar variable z (indicated with Vz ) can be also easily obtained. The clue is to remark that the angular variable θ = ar ctg (y/x), satisfies also the Laplace’s equation, but has an embedded discontinuity at 0 and 2π . Consequently, a fictive slit surface ∂s1 is introduced across the thickness of the membrane and the coordinate Vθ , which runs across the membrane by making a closed loop, is assumed to be discontinuous at the end of the loop itself (Fig. 7). Then, Vθ is chosen to be the solution of the following Laplace’s equation and boundary conditions with [•] standing for the jump of the quantity across the slit:
unc
376
Figure 8 represents the isovalues for both Vθ and Vϕ . Here, we have shown how to obtain either one or three curvilinear coordinates for a general membrane, hollow or not. Usually, as already mentioned, these coordinates cannot be exactly computed, but can be rather easily numerically calculated by the finite element method, so that a good approximation of them is achieved. Although it has not been exploited here, an adaptive finite element mesh refinement could represent a first step to a very interesting pseudo-CAD adaptive approach. This aspect is important because the present method allows assessing the quality of the parametrization in the case of tiny geometrical irregularities with curvature oscillations by simply controlling the approximation error on the Laplace equation itself. In the Appendix we propose a preliminary convergence study.
⎧ ⎨ −div (c∇Vθ ) = 0 inside 0 [Vθ ] = 2π across the cut ∂s1 ⎩ ∂ Vθ ∂n = 0 on ∂e and ∂i
(27)
Again, the above Laplace’s equation can be easily solved by introducing into the finite element formulation a crack-like discontinuity along the ∂s1 slits.
4 Electrical parametrization and mechanics of a thick shell In this section, for illustration purpose how the preceding full system of electric coordinates is used to compute thick shell stresses and strains. If ξs = (ξ1 , ξ2 ) = (Vθ , Vz ) are the tangential coordinates and ξ3 = Vξ3 the thickness variable of the membrane, any point p at its initial configuration can be expressed as similarly as in Eq. (1). When the position is given explicitly from the parameters, as shown in Sect. 2, it is usually easier to get first the
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
TYPESET
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
392 393 394 395 396
397
398 399 400 401 402 403 404 405 406 407 408 409 410 411 412
413 414
415 416 417 418 419 420 421 422
Comput Mech
425
427 428 429 430
431 432 433
434 435
436
437
438
439 440 441
442
443 444 445 446
447
448 449
gp,3
(29)
However here, the introduction of the electric potential provides the inverse mapping from the position to the parameters. Consequently, the contravariant basis is directly obtained first from the gradient of the potentials gsp = ∇ p ξs g3p
(30)
= ∇ p ξ3
Then, with g p = computed as
gp,3
(34)
D p (u) = in ⊗ ∇ p u n
(35)
The weak formulation of the equilibrium of the structure is given by T r [(λT r (ε (u)) I + 2µε (u)) · ε (w)] d V
451
452 453
454
455 456
457
0
∇ p ξ2 ∧ ∇ p ξ3 gp ∇ p ξ3 ∧ ∇ p ξ1 = gp ∇ p ξ1 ∧ ∇ p ξ2 = = n0 gp
(31)
When the shell undergoes a transformation and we write its kinematics, the actual position x, under the classical Mindlin assumptions, becomes x = x 0 (ξs ) + ξ3 R (ξs ) n0 (ξs )
(32)
where x 0 corresponds to the deformed initial position of a point p0 on the middle surface and R is the rotation of the normal fiber n0 . Under small strains hypotheses, the displacement u can be derived u = u0 + ξ3 u1 ∧ n0
=
g1p , g2p , g3p , the covariant basis can be
gp,1 = gp,2
450
where u1s is the tangential component of u1 and Dp (u) stands for the displacement gradient which is computed from
cted
426
∂p ∂ξs ∂p = = n0 ∂ξ3
gp,s =
tensor is given by the following simple expression ε (u) = D p (u) gp,s ⊗s ∇ p ξs + u1s ⊗s ∇ p ξ3
0
f V , w dV
(36)
with f V the volume forces applied to the system, w an admissible test displacement and (λ, µ) the Lame coefficients of the material. Knowing the approximated electric potential coordinates system, the weak form can be easily computed providing the stiffness matrix and the force vector are obtained. For instance, let us consider a membrane with a double curvature such as a portion of the same regular ellipsoid used in this paper. When it is subjected to its own weight, a global deformation occurs as represented in Fig. 9. We have briefly shown here that the electric potential coordinates system can be elegantly used to build the equations and the weak form of shell mechanics. This system directly provides a contravariant coordinates system, which means that the covariant system must be computed from the usual properties of orthogonality from one basis with respect to the other. Obviously, the main advantage of the method is that the coordinates system is numerically obtained via the finite element approximation of the three Laplace’s equations with appropriate boundary conditions. Although not considered here, issues associated with very thin shells and shear
orre
424
covariant basis by simple differentiation
unc
Author Proof
423
pro of
Fig. 8 a The tangential variable Vθ for a curvilinear system of coordinates. b The representation of Vφ in the case of a spherical system of coordinates
(33)
with u0 the middle surface displacement and u 1 the angular rotation vector. Then, it can be shown that the linear strain
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
TYPESET
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
458
459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479
pro of
Comput Mech
Fig. 10 Cross-section of the Drosophila embryo at the onset of gastrulation obtained by confocal 2PEF (two photon excited fluorescence) Fig. 9 A double-curvature shell submitted to its own weight
483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513
5 Application to morphogenesis of Drosophila embryo To conclude this paper, we propose here a particular and rather complex application of the above technique to the morphogenesis of the Drosophila embryo. Morphogenesis is a fundamental process during which an organism achieves its final shape. As many other biological events, it is governed by a perfect synchronization between biology and mechanics. In the last decades, the understanding of its most relevant phases has significantly progressed thanks to interesting numerical and experimental analyses. Furthermore, Drosophila embryo, as other multicellular organisms, presents many similarities with the human embryo, in particular concerning the genetic code. Therefore, a better comprehension of its development may be helpful in providing the necessary insights for the study of tumorogenesis and other diseases [5]. During morphogenesis, and gastrulation in particular, several morphogenetic movements take place which lead to significant changes in shape. At this moment of the development, the form of the embryo is very similar to the one of a bean (Fig. 10) with a dorsal region flatter than the ventral one and two rounded poles, the anterior pole (AP) and the posterior pole (PP). The embryonic tissue, the blastoderm, is composed by a single layer of about 6,000 epithelial columnar cells that are arranged in a series of almost circular arrays along the horizontal axis. Each cell has the apical surface in contact with the external semi-rigid vitelline membrane, which surrounds the entire structure. The basal surface instead faces the yolk, a viscous and slightly compressible fluid filling the internal region of the embryo. The apico-basal axis of each cell is aligned with the axis of radial symmetry of the embryo. The thickness of the embryo is then in the order of some microns which gives to it the aspect of a
cted
482
orre
481
or membrane locking can be considered as in the classical framework.
unc
Author Proof
480
three-dimensional membrane (Fig. 10). For this reason, we think that the new technique previously introduced can be profusely employed to this specific case in order to reproduce the behavior of the system. Furthermore, the elementary deformations that trigger the morphogenetic movements have been identified through experimental observations by biologists, so that they can be easily plugged into the analytical formulation. In our previous work, we proposed a finite element model to simulate multiple morphogenetic movements of a simplified ellipsoidal Drosophila embryo [3]. Even though the results are very interesting and consistent compared to the reality, the theoretical approach requires several mathematical manipulations that may become very difficult when the shape of the embryo deforms and becomes more complex. Therefore, the main advantage of the new method lies in a more intuitive and reduced analytical formulation, which allows describing much more general three-dimensional membranes. 5.1 Morphogenetic movements computed from the electric coordinates system of the ellipsoid During morphogenesis, each cell tends to deform itself as if it was completely free, what is called the active deformation. But actually, the cells are constrained by external or internal components of the system (i.e. vitelline membrane and yolk) and they are also in contact with each other. Therefore, the active strains may generate some incompatibilities with respect to the continuity of the system, such as a superposition of the cellular domains. Then, the cells react and rearrange themselves by an elastic passive response in order to guarantee the boundary conditions imposed by the biological system. From a mechanical standpoint, it has been so far very difficult to distinguish between active and passive deformations. This is mainly due to the fact that embryonic cells do not act alone, but they are strongly coordinated with their neighbors. Nevertheless, from experimental observations, biologists have been able to identify some elementary
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
TYPESET
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532
533 534
535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550
Comput Mech
554 555 556 557 558 559 560 561 562 563 564 565 566
568 569 570 571 572 573 574 575 576 577
578
579 580 581 582 583
584
When the elementary active deformations take place, any point p moves to an intermediate configuration x, different for each morphogenetic movement. Deriving x will yield to the expression of the active deformation gradient F a , which, according to Eq. (5), is defined as follows [22] Fa =
∂x ∂x ∂x ⊗ ∇ p Vθ + ⊗ ∇ p Vz + ⊗ ∇ p Vξ3 ∂ Vθ ∂ Vz ∂ Vξ3
(38)
585
586 587 588
589
590 591 592 593 594 595 596
597
As mentioned previously, the total deformation gradient F, the composition of the active and the passive deformations, is determined as F = Fm Fa
unc
Author Proof
567
where f s indicates the surface forces and w is any displacement test function. Even though they are not geometrically modeled here, the vitelline membrane and the internal yolk are taken into account and implemented in the right hand side term of Eq. (40), as similarly as in our previous work [3], by contact and constant volume boundary conditions, respectively. This choice provides a drastic decrease of the number of d.o.f. and a largely reduced computational cost. Here, the electric potential method will only serve for the evaluation of the active gradient F a , whereas the blastoderm will be modeled as a three-dimensional continuum domain.
pro of
553
strains that occur during the biological movements and lead to the final specific configurations at different stages of morphogenesis [9,15,17–19,23,25]. In order to take into account this aspect, we use here the deformation gradient decomposition [22], a method very often introduced in biomechanics and already proposed in our previous work [3] and by other authors [4,8]. This approach allows splitting the final deformation of an infinitesimal material element of the structure into two parts: a non-stress generating deformation Fa (the active deformation) and a stress-generating deformation Fm (the passive deformation). The Drosophila embryo is approximated by a regular ellipsoid with real dimensions. Thus, it is 500 µm in length along the anterior–posterior axis and 150 µm along the two transversal directions. The thickness of the cellular layer has been chosen uniform and equal to 15 µm. In our preceding publication [3], the ellipsoid as well as the covariant and the contravariant bases, the middle surface and the curvature have been analytically computed. Here instead, these quantities are obtained by the electric potential technique presented in Sect. 3. Using the complete curvilinear basis found in Sect. 3.2, p0 (Vθ , Vz ) [Eq. (10)] is a point on the middle surface of the ellipsoid such that Vξ3 (x, y, z) = 0, with Vθ and Vz the associated coordinates to V 1 and V 2 , respectively. We can define any point p Vθ , Vz , Vξ3 through the thickness of the ellipsoid as p Vθ , Vz , Vξ3 = p0 (Vθ , Vz ) + Vξ3 n0 (Vθ , Vz ) (37)
5.1.1 Ventral furrow invagination (VFI)
(39)
To describe the isotropic and elastic behavior of the embryonic tissue, we use a simple Saint-Venant model [22] and we choose the Young’s modulus E = 100 Pa [26] and the Poisson’s ratio ν = 0.45. The weak form of the equilibrium condition in the initial configuration 0 is expressed through the first Piola Kirchhoff stress π as follows T (40) w, f s J F−T (n) d S T r π Dp w d V =
m
∂m
The VFI is one of the most studied among the morphogenetic movements. Several numerical models, in two or three dimensions, have been proposed in literature in the last decades [1,4,7,8,14,21,24]. It is triggered by a series of changes in shape of single cells in the areas of active movement [17,19,23]. The process has been separated into several successive steps: flattening of the apical surfaces at the ventral region of the cellular blastoderm [17,23]; constriction of the apical surface of some spread cells [1]; lengthening along the apico-basal axis of the cells within the furrow [9]; shortening back of the cells to the original length once the ventral furrow is formed, while keeping the apical ends constricted [9,18]; lateral pushing exerted by lateral and dorsal cells, which facilitates and reinforces the invagination of the ventral furrow [4] (Fig. 11a–c). Practically, the cells change their shape from a cylindrical to a more trapezoidal form, inducing the curvature necessary for the invagination of the blastoderm into the yolk. Biologists seem to believe that apical constriction and apico-basal elongation are the main responsible deformations for VFI, while the other three deformations are an elastic response of the cells. Partly, we have been able to confirm this hypothesis with our previous work [3]. We observed in fact, that introducing only apical constriction as active deformation, the invagination of the cellular layer still occurs. Therefore, for the present work, we still keep this strain as the only initiator of the movement. The deformation is applied on a limited region of the embryo, determined by a regularized level set function, which corresponds to the actual active area where the genetic process takes place. Therefore, the intermediate position xVFI Vθ , Vz ,Vξ3 for this specific case, where the variables Vθ , Vz , Vξ3 have been computed as described in Sect. 3.2, can be written as xVFI = p0 V θ , Vz + Vξ3 n0 V θ , Vz (41) with
V θ = Vθ + α ac
TYPESET
600 601 602 603 604 605 606 607 608 609
611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643
644
645
2Vζ3 m (Vθ ) h
(42)
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
598 599
610
cted
552
orre
551
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
646
Comput Mech
5.1.2 Germ band extension (GBE)
660
Author Proof
pro of
Once gastrulation has been completed, the germ band, located at the ventral region of the embryo, starts to extend. This movement progresses to bring the posterior tip of the germ band to about 75% egg length (0% length coincides with the posterior pole) (Fig. 11d–f). The elongation of the tissues is due to a convergence–extension movement of a population of cells at the central–lateral region of the embryo. This process is triggered by an intercalation of the cells that interpose themselves between their dorsal or ventral neighbors, resulting in a decrease of the number of cells along the dorsal–ventral axis and in an increase of the number of cells along the anterior–posterior axis [15]. Some interesting models have been proposed in the literature in the last decades [25]. As for the previous work, we do not precisely simulate the intercalation since the cells are not geometrically modeled. However, we propose a uniform movement of compression–extension tangential to the middle surface of the blastoderm. As for the VFI, the active deformations are introduced on a limited region of the embryo, in this case corresponding to the germ band at ventral region. The interGBE Vθ , Vz , Vξ3 is equal this time to mediate position x xGBE = p0 V θ , V z + Vξ3 n0 V θ , V z (43)
cted
Fig. 11 a–c Experimental images of cross-sections of the Drosophila embryo during VFI. d–f Sketch of the successive steps of GBE
where
650 651 652 653 654 655 656 657 658 659
Fig. 12 a–c Successive steps of the numerical simulation for the VFI using the complete formulation presented in Sect. 3.2
V θ = 1 + αθGB E Vθ V z = 1 + αzGB E Vz
DISK
LE
668 669 670 671 672 673 674 675 676 677 678 679 680 681
682
(45)
685
with αθGBE and αzGBE the amplitudes, respectively, of the shortening from the dorsal to the ventral region of the embryo and of the lengthening along the anterior–posterior axis. The results for this simulation are very similar to the ones previously obtained [3] (Fig. 13). Actually, it is possible to clearly observe the convergence of the embryonic tissue from the dorsal to the ventral region of the embryo together with the extension towards both the AP and the PP. Additionally, the characteristic vortex movement, particularly at the PP,
TYPESET
666 667
684
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
664 665
(44)
orre
649
where is the amplitude of the apical constriction. The periodic function m (Vθ ) mimics the material cells by which the active region is subdivided and it is the same used in our previous work [3]. In Fig. 12, the phases of the numerical simulation are shown. The amplitude of the invagination is somewhat smaller than the one obtained in our previous analysis [3]. However, it has to be said that the main objective of the present study is to propose an alternate method to describe the biological event. Thus, the evaluation of the results is so far purely qualitative. Nevertheless, we think that the application on a rather complex biological geometry such as the one of the Drosophila embryo constitutes a helpful step to validate the approach.
663
683
unc
647 648
α ac
661 662
CP Disp.:2010/11/4 Pages: 16 Layout: Large
686 687 688 689 690 691 692 693 694
Comput Mech
698 699
Author Proof
700
701 702
703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734
Fig. 14 The new geometry of the Drosophila embryo obtained by simply scaling the original shape of the regular ellipsoid
5.2 A more accurate representation of the embryo initial geometry In the previous section, as well as for our preceding work [3], the ellipsoid has been employed to model the Drosophila embryo and it is actually a good representation of the real geometry. Besides, as other regular three-dimensional forms (spheres, cylinders, etc.), when one of the dimensions is sufficiently small it can be easily parametrized, as shown in Sect. 2. The analytical description becomes more complicated when the geometries are more complex and present irregularities such as concavities or convexities. This could be the case of the Drosophila embryo, for which the ellipsoid is a convenient but not fully accurate approximation. Actually, some restrictions may be noticed. For instance, the ellipsoid has a symmetric shape with respect to the two principal axes. Nevertheless, in the Drosophila embryo (Fig. 10) the PP is more rounded than the AP and the dorsal region is flatter (almost concave) than the ventral region, which is however more rounded, which implies that the form of the embryo is rather similar to a bean. Taking into account these two aspects may be relevant in order to obtain more consistent results. First, the rounding up of the PP would be fundamental during GBE, while the lengthening of the ventral tissue takes place. A sharpened pole, as in the case of the ellipsoid, may instead inhibit this movement. Furthermore, the GBE towards the posterior pole leads to a series of little invaginations at the dorsal region of the embryo, which is probably helped by the flattening of this area. Therefore, in the case of the ellipsoid, this phenomenon could be hindered twice. In this section, we consider a modified ellipsoid, which has been obtained by simply scaling the original one. As we can observe in Fig. 14, the new shape presents a more rounded PP and a flatter dorsal region. Its dimensions are slightly larger
than the ellipsoid previously used. The major axis is 500 µm long, while the cross axes CE and DF are, respectively, 175 and 165 µm long. As for the case of the ellipsoid, we propose the simulation of two morphogenetic movements: the VFI and the GBE. In Sects. 5.1.1 and 5.1.2 we have shown how the elementary deformations are directly introduced into the analytical formulation through the complete curvilinear basis obtained in Sect. 3.2. Here, we have decided to use the reduced formulation presented in Sect. 3.1 to simplify the problem. Therefore, we have calculated the electric potential (Fig. 15a), which allows us to determine the three vectors V 1 , V 2 , V ξ3 as shown in Fig. 15b, c and d, respectively. As mentioned above, this approach is valid as long as Fn,m is given a priori and the elementary deformations can be simply described. This is the case here. Actually, the VFI can be represented by a bending in the direction of V 1 , while the GBE corresponds to a compression–extension movement in the directions of V 1 and V 2 , respectively. Therefore, the active deformation gradients FaVFI and FaGBE , respectively, for the VFI and the GBE read
cted
697
orre
696
can be noticed. Once more, we can conclude that our new approach can be considered as an efficient tool to describe specific biological events. Furthermore, as mentioned in the previous section, here we only show the versatility of the new technique from a qualitative point of view. More quantitative results will be detailed in a forthcoming paper.
unc
695
pro of
Fig. 13 Results for the numerical simulation of the GBE; the complete covariant basis has been employed as shown in Sect. 3.2
FaVFI = I + α ac Vξ3 V 1 ⊗ V 1 FaGBE = I + αθGBE V 1 ⊗ V 1 + αzGBE V 2 ⊗ V 2
DISK
LE
738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755
756
(47)
757
Some differences may be noticed, in particular for the VFI where the amplitude of the invagination is again smaller than the one obtained in our previous work (Fig. 16). Nevertheless, the results are still consistent from a qualitative point of view, especially for the simulation of the GBE (Figs. 17, 18). The main goal of this section of the work is to show that, even with a partial parametrization, we are able to describe irregular geometries and provide nice simulations of such
TYPESET
736 737
(46)
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
735
CP Disp.:2010/11/4 Pages: 16 Layout: Large
758 759 760 761 762 763 764 765
Comput Mech
Author Proof
pro of
Fig. 15 a Computation of the electric potential Vξ3 . b The normal vector n0 . c, d The tangential vectors V θ and V z , respectively
orre
cted
Fig. 16 Numerical simulation of the VFI for the new geometry of the Drosophila embryo
767
complex biological events. The detailed influence of the new geometry will be presented also in a forthcoming paper.
768
6 Conclusions and perspectives
769 770 771 772 773 774 775 776 777 778 779 780
In this paper, we propose an alternative and original method to analytically describe thick biological membranes with complex geometries. The approach is based on the strong connection between the electric potential and the domain inside which it is calculated. The main idea is to parametrize the so-called normal coordinate of any general thick shell using the electric potential and to define then an associated curvilinear basis. Two forms of the formulation are proposed: a reduced form and a complete form, which incorporates the tangential variables also. In order to illustrate our approach, we propose here an interesting application to the Drosophila embryo, whose geometry may become very complex during
unc
766
Fig. 17 Results for the numerical simulation of the GBE for the new scaled ellipsoid
morphogenesis. The individual simulations of two morphogenetic movements (ventral furrow invagination and germ band extension) are presented together with relative consistent results.
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
TYPESET
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
781 782 783 784
Comput Mech
pro of 793
Appendix: A preliminary convergence analysis
787 788 789 790 791
794 795 796 797 798 799 800 801
802
803 804 805 806
orre
786
cted
792
Finally, the present technique may be very useful when combining together several morphogenetic movements into a concurrent simulation since it allows taking into account peculiar aspects so far not considered as it will be presented in a further publication. Furthermore, coupling electrics and mechanics would allow later to take into account the potential electrical signals developing inside the embryonic cells, like already considered in heart modeling.
785
Fig. 19 Two-dimensional annular domain
In this appendix we consider a first convergence analysis on a simple two-dimensional model. An extended study will be provided in a future publication. We consider an annular domain , as represented in Fig. 19, with radius r such that ri ≤ r ≤ re , with ri and re the internal and the external radius, respectively. We look for an approximate solution of the Laplace’s equation with unknown u ⎧ ⎪ ⎨ u = 0 ∂u ∂n e = σe ⎪ ⎩ ∂u = σi ∂n i
on on ∂e on ∂i
unc
Author Proof
Fig. 18 Example of completely deformed ellipsoid. Specifically, this is the shape of the Drosophila embryo once the VFI has occurred. With the new technique we are able to compute the electric potential Vξ3 (a) and to obtain the three orthogonal vectors n0 (b), Vθ (c) and Vz (d)
The compatibility condition requires that the fluxes satisfy σi ri + σe re = 0
(A2)
The exact solution depends on r only and it is given by
(A1)
where n e and n i are the external and the internal outwards normal vectors, σe and σi the external and the internal fluxes, ∂e and ∂i the external and the internal boundaries of the two-dimensional domain (Fig. 19).
u = a ln(r ) + b
(A3)
with b an undetermined constant and a = σe re . Here, the natural coordinates are of course the usual polar variables r and θ , but we will apply the electrical method to get an approximation of them and of the solution of the Laplace’s equation itself.
TYPESET
DISK
LE
808
809
810
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
807
CP Disp.:2010/11/4 Pages: 16 Layout: Large
811
812 813 814 815 816
Comput Mech
pro of
ments and the electrical approach, we find ah = −0.899402 and for the refined mesh ah = −0.899168, which means that a very good approximation has already been obtained even with the crude mesh. In the case of a solution with θ dependancy of the flux, e.g. σi = sin θ , similar good results have been obtained. Therefore, although it is a first attempt to qualify the approximation properties of the method, it is believed that from this standpoint the results are very encouraging. Obviously, further extended numerical analysis is planned for more complex examples such as membranes and shells.
Fig. 20 c, d The isovalues of rh and θ for a coarse (a) and a refined (b) mesh
819
820 821 822 823 824 825 826 827 828
829
830 831 832
833
If we choose ri = 0.9 and re = 1, the isovalues of rh and θh are given for a coarse mesh (Fig. 20a, c) and a refined mesh (Fig. 20b, d) with quadratic elements. We see that rh raises correctly almost linearly from 0.9 to 1 and θh from 0 to 2π even for a crude mesh. Let us now evaluate the approximate solution u h of u given by Eq. (A3). In the spirit of the electrical method and its application to shells or membranes, we simply assume a trial function given by u h = ah ln(rh ) + bh
834
836 837 838 839 840 841
(A5)
where bh is again an undetermined constant. We compute ah by a Galerkin approach, which means that it is the solution of the underneath variational equation σe uˆ h dl + ah ∇u h , ∇ uˆ h d S = σi uˆ h dl (A6)
835
1. Alberch O et al (1981) The mechanical basis of morphogenesis. Epithelial folding and invagination. Dev Biol 75:446–462 2. Allain JM (2005) Instabilité des membranes lipidiques inhomogènes. Implications biologiques. PhD thesis, Université Paris VII 3. Allena R, Mouronval A-S, Aubry D (2010) Simulation of multiple morphogenetic movements in the Drosophila embryo by a single 3D finite element model. JMBBM 3:313–323 4. Barret K et al (2007) A deformation gradient decomposition method for the analysis of the mechanics of morphogenesis. J Biomech 40:1372–1380 5. Bate M, Arias AM (1993) The development of Drosophila melanogaster. Cold Spring Harbor Laboratory Press, Cold Spring Harbor 6. Brown ML (1983) Scalar potential in multiply connected beam. Int J Numer Methods Eng 20:665–680 7. Clausi DA, Brodland GW (1994) Embryonic tissue morphogenesis modeled by FEM. J Biomech Eng 116:146–156 8. Conte V et al (2007) A 3D finite element model of ventral furrow invagination in the Drosophila melanogaster embryo. J Mech Behav Biomed Mater 1:188–198 9. Costa M et al (1993) Gastrulation in Drosophila: cellular mechanisms of morphogenetic movements. In: Bate M, Martinez-Arias A (eds) The development of Drosophila. Cold Spring Harbor Laboratory Press, New York, pp 425–466 10. Courant R, Hilbert D (2007) Methods of mathematical physics, vol 2. WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 11. Green AE, Zerna W (1968) Theoretical elasticity. Oxford University Press, Oxford 12. Greenbaum A, Greengard L, McFadden GB (1993) Laplace’s equation and Dirichlet–Neumann map in multiply connected domains. J Comput Phys 105:267–278 13. Hunter PJ, Nielse PM, Smaill BH, LeGrice IJ, Hunter IW (1992) An anatomical heart model with applications to myocardial activation and ventricular mechanics. Crit Rev Biomed Eng 20(5–6):403–426 14. Jacobson AG et al (1986) Neurulation and cortical tractor model for epithelial folding. J Embryol Exp Morph 96:19–49 15. Keller R et al (2000) Mechanisms of convergence and extension by cell intercalation. Philos Trans R Soc Lond B Biol Sci 355(1399):897–922 16. Kellog OD (1953) Foundations of potential theory. Berlin Verlag Von Julius Springer, Dover 17. Leptin M, Grunewald B (1990) Cell shape changes during gastrulation in Drosophila. Development 110:73–84 18. Leptin M (1995) Drosophila gastrulation: from pattern formation to morphogenesis. Annu Rev Cell Dev 11:189–212
cted
818
Let rh and θh be such electrical coordinates which consequently satisfy the following equations ⎧ on ⎧ ⎪ ⎪ θh = 0 ⎪ ⎨ ∂θh = 0 ⎨ rh = 0 on on ∂e ∂n e rh = re on ∂e (A4) ∂θh ⎪ ∂n = 0 ⎩ on ∂i i ⎪ r h = ri on ∂i ⎪ ⎩ [θh ] = 2π across ∂s
∂e
∂i
with the test function uˆ h = aˆ h ln(rh ) + bˆh . For the numerical application we choose σi = 1 and σe = −0.9. The exact value for a, for the polar coordinates case, is thus σe re = −0.9, which is also the value provided for ah in Eq. (A5) with rh replaced by r . For the trial function as expressed in Eq. (A5), ah is equal to −0.898337 which is rather close to the exact solution. When u h is computed by the finite element method, for the coarse mesh with quadratic ele-
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
TYPESET
DISK
LE
843 844 845 846 847 848 849 850 851 852
853
orre
817
unc
Author Proof
References
842
CP Disp.:2010/11/4 Pages: 16 Layout: Large
854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899
Comput Mech
901 902 903 904 905 906 907 908
24. Taber LA (2004) Non linear elasticity: applications in biomechanics. World Scientific Publishing, Singapore 25. Weliky M, Oster G (1990) The mechanical basis of cell rearrangement. Development 106:372–386 26. Wiebe C, Brodland GW (2005) Tensile properties of embryonic epithelia measured using a novel instrument. J Biomech 38: 2087–2094 27. Zienkiewicz OC, Taylor RL (1987) The finite element method. McGraw-Hill, London
unc
orre
cted
Author Proof
909
19. Leptin M (1999) Gastrulation in Drosophila: the logic and the cellular mechanisms. EMBO J 18:3187–3192 20. Niordsson FI (1985) Shell theory. North-Holland, Amsterdam 21. Pouille PA, Farge E (2007) Hydrodynamic simulation of multicellular embryo invagination. Phys Biol 5:015005 22. Smith DR (1993) An introduction to continuum mechanics. Kluwer, The Netherlands 23. Sweeton D et al (1991) Gastrulation in Drosophila: the formation of the ventral furrow and posterior midgut invaginations. Development 112:775–789
pro of
900
123 Journal: 466 Article No.: 0551 MS Code: CM-10-0284
TYPESET
DISK
LE
CP Disp.:2010/11/4 Pages: 16 Layout: Large
910 911 912 913 914 915 916 917 918
