1. Introduction

Extracellular, diffusible signals provide critical pieces of information to cells in biological systems. Often, the concentration of these signals conveys additional, spatial information. For example, during animal development, morphogens prompt concentration-dependent cell fate decisions to pattern tissues ^[1,2,3], and chemoattractants provide directional information to migratory cells ^[4,5,6,7,8]. Work in vitro clearly demonstrates these principles; however, the situation in vivo is more complex. While there is debate about how diffusible protein-based signals traverse through tissues to create gradients (see examples in ^[1]), one simple possibility is that such molecules move between and around cells. In this case, the signals are likely to be hindered by cell packing and to pool in interstitial spaces, which could effectively change their perceptible concentration and gradients. Even in cell-dense animal tissues, small spaces or gaps can exist between non-epithelial cells, through which signals can move. Such domains have been characterized in amphibian embryo ^[9,10], and in mammalian brain, where extracellular spaces may make up 20% of the overall tissue ^[11]. Although extracellular matrix components may occupy some of this interstitial space and detain diffusible molecules, often not all space is occupied, leaving significant areas unoccupied/fluid-filled for diffusion.

To examine how the interstitial space in tissue affects signal distribution, we modeled the elegant case of the Drosophila melanogaster egg chamber. Fly oogenesis has been extensively characterized and has provided a tractable model to elucidate the genetic regulation of coordinated cell development, tissue patterning, stem cell maintenance, and cell migration ^{[12,13,14,15]}. In Drosophila, eggs develop within an egg chamber, which is comprised of an oocyte, fifteen large, germline sister cells called "nurse cells" that foster oocyte growth, and an outer, somatic, epithelial cell layer that is essential for patterning and eggshell formation ^[16]. Although nurse cells are space-filling, their large size and constraints created by cell junctions result in notable spaces between these cells and between them and the enveloping follicle cells. Among the epithelial cells, a set of migratory cells—the border cells—arise during mid-oogenesis, and migrate as a cluster between nurse cells in response to extracellular chemoattractant cues at oogenesis stages 9 and 10 ^[17,18] (Figure 1).

Motile cells often move directionally in response to gradients of chemoattractant or chemorepellant molecules. Border cell migration to the oocyte requires several chemoattractants, which bind and activate two conserved receptor tyrosine kinases (RTK) on the border cells: Epithelial Growth Factor Receptor (EGFR) and a receptor equally related to mammalian Platelet Derived Growth Factor Receptor and Vascular Endothelial Growth Factor Receptor, which is called PVR (reviewed in ^[17,18]). Both receptors appear to be redundantly required for directed migration to the oocyte at stage 9 ^{[19,20,21,22]}, and EGFR is additionally required for migration to the dorsal side at a later stage ^[23]. Multiple ligands exist for each receptor and may act redundantly as guidance cues ^[19,20,24]. These factors are produced at the oocyte, but may also be generated from some somatic epithelial cells ^[20,24,25]. Since the RTK ligands function as chemoattractants, the simplest model is that these cues are released from the oocyte surface and evolve into a gradient across the egg chamber, driving clustered cell migration. It has been shown experimentally that uniform, local overexpression of chemoattractant or reduced RTK activity disrupts migration ^{[19,20,21,22,23,24]}, and ectopic chemoattractant expression can mis-direct the cells ^[20,24], demonstrating that the concentration and gradient are both important. However, chemoattractants are challenging to detect, observe, or measure in vivo, so their precise distribution in time and space is unknown. Moreover, we postulated that the gradient may be altered by hindrance due to the large nurse cells that are packed into this domain. Several questions surround the characteristics of chemoattractant in the extracellular regions within the egg chamber:

(1) How does the heterogeneous egg chamber environment impact the local quantitative values of chemoattractant?

(2) How does secretion location impact chemoattractant localization?

(3) How does the extracellular space between the nurse cells impact the gradient of chemoattractant in a realistic geometry?

(4) How does the gradient in this environment compare without the nurse cell heterogeneity?

(5) How does migration timing differ due to the chemoattractant distribution in this geometry?

To address how cell packing in a tissue influences signaling, we generated a geometrically-realistic model of the egg chamber then mapped the distribution of chemoattractants in this domain using a reaction-diffusion system. Using a biophysically-relevant range of parameters, our model indicates that guidance cues diffusing only in interstitial space create a more dramatic concentration gradient than if they moved throughout the whole egg chamber unimpeded by nurse cells. In addition, the gradient is established much more quickly when large cells are packed in the domain, which may be important to allow development to proceed on time. Interestingly, the directionality and the magnitudes of gradient vectors along the migratory route varied by position in the egg chamber, which may explain some aspects of the cell migration observed in vivo. Thus, our model reveals several unanticipated features of diffusible signaling within complex animal tissue and provides a realistic input that can inform models of cell migration in response to chemoattractants ^[26,27,28].

---

2. Materials and methods

---

2.1. Creating the extracellular domain within the egg chamber

In previous work, we modeled border and polar cell cluster migration in the egg chamber of Drosophila melanogaster ^[26]. The Stonko model represents each cell (epithelial follicle, border, polar, nurse cells) in the anterior egg chamber at stage 9 with one or more Identical Math Cells (IMCs) depending on each cell's relative size. Several forces impact the motion of the cluster, and one of these forces is chemotactic, the response to chemoattractant. In our previous work, the gradient of chemoattractant was assumed to be constant and uniform from low at the anterior end to high at the oocyte, and the border cells moved steadily up the gradient.

Here, we model the diffusion and reaction of chemoattractant in the extracellular domain within an egg chamber using a linear parabolic partial differential equation. Chemoattractant is secreted from the oocyte boundary, and in some simulations from other areas of the epithelia, as flux boundary conditions. Degradation of the chemoattractant is assumed within the extracellular space, while some uptake is assumed on the surface of the nurse cells and taken as a mixed boundary condition.

To obtain the geometry of the extracellular domain, we imported the IMCs from the force-based Stonko model, and generated the convex hull of the center of the IMCs representing the initial boundary of the anterior half of an egg chamber (the border, polar and epithelial cells, as well as the oocyte). Then for each nurse cell, we generated an alpha shape containing the centers of the IMCs, and used the volume formed by the outer edges of these boundaries to exclude nurse cells from the extracellular domain (Figure 2). In this model, the egg chamber is approximately a paraboloid of length 409.5 μm ending in the oocyte, a circle of radius 146.25 μm, and centered around y = z = 146.25 μm. This is an upper bound for a Drosophila stage 9 egg chamber size. We also consider a domain with length and radius half that size for comparison.

---

2.2. Model of secretion and diffusion in the egg chamber extracellular space

Chemoattractants originate from a source at the oocyte, and potentially from the epithelium, and diffuse within the egg chamber, increasing the concentration at the anterior across the egg chamber. The polar and border cells exit from the anterior epithelium and respond to chemoattractants by migration towards the oocyte ^{[19,20,23,24]}. Included below are the equations used in our model of chemoattractant diffusion and reaction.

$\frac{{\partial {\rm{u}}}}{{\partial {\rm{t}}}}{\rm{ = D}}\left( {\frac{{{\partial ^{\rm{2}}}{\rm{u}}}}{{\partial {{\rm{x}}^{\rm{2}}}}}{\rm{ + }}\frac{{{\partial ^{\rm{2}}}{\rm{u}}}}{{\partial {{\rm{y}}^{\rm{2}}}}}{\rm{ + }}\frac{{{\partial ^{\rm{2}}}{\rm{u}}}}{{\partial {{\rm{z}}^{\rm{2}}}}}} \right){\rm{ - ku}}$

(1)

${\left. {{\rm{ - D}}\frac{{\partial {\rm{u}}}}{{\partial {\rm{n}}}}} \right|_{{\rm{boundary}}\;{\rm{source}}\;{\rm{of}}\;{\rm{chemoattractant}}}}{\rm{ = \sigma }}$

(2)

${\left. {{\rm{ - D}}\frac{{\partial {\rm{u}}}}{{\partial {\rm{n}}}}} \right|_{{\rm{boundary}}\;{\rm{of}}\;{\rm{nurse}}\;{\rm{cells}}}}{\rm{ = - }}\phi {\rm{u}}$

(3)

${\left. {{\rm{ - D}}\frac{{\partial {\rm{u}}}}{{\partial {\rm{n}}}}} \right|_{{\rm{all}}\;{\rm{other}}\;{\rm{boundaries}}}}{\rm{ = 0}}$

(4)

The base of our model, Eq 1, is the diffusion-reaction equation for the chemoattractant concentration, u, where D is the diffusion coefficient, and k (1/s) represents the rate at which the chemoattractant breaks down in the extracellular space of the egg chamber. The boundary conditions, Eqs 2–4, model the behavior of diffusion at different types of surfaces within the egg chamber. These are Neumann conditions, modeling the derivative of the concentration in the direction normal to the boundary, in this case, the flux. Eq 2 represents the secretion of chemoattractants with parameter σ (pM μm/s) being the value of the source, taken to be constant, while Eq 3 models the uptake of the chemoattractant, which depends on the concentration at the boundary, with parameter ϕ (μm/s) being the rate constant. The no-flux boundary condition in Eq 4 indicates that the boundaries at all other points do not let any chemoattractant in or out. The boundary conditions can easily be changed to model secretion or uptake at any boundaries to test different hypotheses.

Typically we simulated for 18000 seconds (5 hours), as this timespan corresponds with the developmental time line of the oogenesis stages during which border cells form and migrate in D. melanogaster. We varied several of the critical model parameters from their default values: D = 1 μm² s^–1, k = 1 × 10^–4 s^–1, ϕ = 0 μm s^–1, and σ = 100 pM μm s^{–1 [29]}. In our tests, unless otherwise stated, the chemoattractant is secreted from the region of the epithelium where the x-value is greater than 300 μm. The chemoattractant is also secreted from the face of the oocyte, located at x-value 409 μm. Given that these default choices can only be estimates, we varied the parameter values to check the sensitivity of our measures after nondimensionalizing the system.

We non-dimensionalize to reduce the number of independent variables. Let x = LX, y = LY, z = LZ for some system length and non-dimensional spatial variables, X, Y and Z. Let t = aτ for some system time scale a and non-dimensional time variable. Let u = u₀U for some system concentration scale, u₀, and non-dimensional concentration U. Applying these changes of variables, Eqs 1–4 become:

$\frac{{\partial U}}{{\partial \tau }} = \frac{{Da}}{{{L^2}}}\left( {\frac{{{\partial ^2}U}}{{\partial {X^2}}} + \frac{{{\partial ^2}U}}{{\partial {Y^2}}} + \frac{{{\partial ^2}U}}{{\partial {Z^2}}}} \right) - akU$

($1a$)

${\left. { - \frac{{Da}}{{{L^2}}}\frac{{\partial U}}{{\partial n}}} \right|_{boundary\;source\;of\;chemoattractant}} = \frac{{a\sigma }}{{{u_0}L}}$

($2a$)

${\left. { - \frac{{Da}}{{{L^2}}}\frac{{\partial U}}{{\partial n}}} \right|_{boundary\;of\;nurse\;cells}} = - \frac{{a\phi }}{L}U$

($3a$)

${\left. { - \frac{{Da}}{{{L^2}}}\frac{{\partial U}}{{\partial n}}} \right|_{all\;other\;boundaries}} = 0$

($4a$)

We pick L = 400 μm, the length of the egg chamber, to fix the domain. We scale time by the domain decay constant a = 1/k. We let u₀ = 1 pM so the threshold to register chemoattractant, u = 1 pM is fixed. With these changes the non-dimensional equations become:

$\frac{{\partial U}}{{\partial \tau }} = {D_{eff}}\left( {\frac{{{\partial ^2}U}}{{\partial {X^2}}} + \frac{{{\partial ^2}U}}{{\partial {Y^2}}} + \frac{{{\partial ^2}U}}{{\partial {Z^2}}}} \right) - U$

($1b$)

${\left. { - {D_{eff}}\frac{{\partial U}}{{\partial n}}} \right|_{boundary\;source\;of\;chemoattractant}} = {\sigma _{eff}}$

($2b$)

${\left. { - {D_{eff}}\frac{{\partial U}}{{\partial n}}} \right|_{boundary\;of\;nurse\;cells}} = - {\phi _{eff}}U$

($3b$)

${\left. { - {D_{eff}}\frac{{\partial U}}{{\partial n}}} \right|_{all\;other\;boundaries}} = 0$

($4b$)

Where ${D_{eff}} = D/(k \cdot {L^2})$, ${\sigma _{eff}} = \sigma /\left({1pM \cdot k \cdot L} \right)$, and ${\phi _{eff}} = \phi /\left({k \cdot L} \right)$. The typical simulation time will then be about = (18000 s) (1 × 10^–4 s^–1) = 1.8.

To get a preliminary idea of how this chemoattractant distribution affects clustered cell migration, we created simulations based on our previous work ^[26], but now allow for a non-constant gradient. So the gradient [U_x, U_y, U_z]^T experienced by each border cell may either be a constant [1, 0, 0]^T or interpolated from the solution to Eqs 1–4.

---

2.3. Numerical implementation of model

Due to the complex shape and boundary conditions, we used the Matlab partial differential equation toolbox and function pdesolve to solve the system 1–4 using finite elements. This solver implements the finite element method to approximate a solution. The maximal element size input parameter was hmax = 100/L. Lower values of hmax did not appreciably alter the output values. Simulations were done in Matlab Version 2016a. Some simulations were performed on the HPCF at UMBC (www.hpcf.umbc.edu).

---

3. Results

---

3.1. Chemoattractant can spread across the egg chamber in minutes

It is understood that a chemoattractant is necessary for the directed migration of the polar and border cell cluster (reviewed in ^[17,18,30]). Therefore it is important, when considering diffusion, to know when the chemoattractant reaches a critical level at the anterior end of the egg chamber, where these cells initially reside. For this reason, we chose the farthest point from the source of secretion to test for a substantial concentration of chemoattractant. We chose the point (20 μm/L, 146.25 μm/L, 146.25 μm/L), located on the anterior end of our egg chamber central in the y and z directions. Before simulating how long it would take for the chemoattractant to reach the anterior side of the egg chamber, we defined a minimal concentration of chemoattractant that anterior cells might detect and that could trigger directed migration. For this we chose 1 pM as a baseline for considering what has sufficient chemoattractant concentration. We show the time to reach this threshold in Figure 3 with nondimensional parameters but also include values for dimensional parameters for comparison (these results are listed in Table 1, see Supplementary movie 1). Further results will appear in dimensional values.

To determine reasonable value ranges for dimensional parameters, we varied them one at a time while holding the others fixed. We first tested the time needed for the chemoattractant to be perceived by the motile cells ("time across") using various D values ranging from 0.01 μm² s^–1 to 100 μm² s^–1 on a logarithmic scale, as listed in Table 1. The "time across" is biologically important since if the migratory cells are not triggered at the right developmental stage, they will not be able to reach their required destination. With the values of 0.01 μm² s^–1 and 0.1 μm² s^–1, the diffusion of the chemoattractant did not reach the end of the egg chamber by five hours, roughly the time to transition between stages 8 and 9 (when the cells are specified then must migrate), leading us to believe that D needs to be greater than 0.1 μm² s⁻¹ to elicit the response observed in vivo. This supports a D value at least on the order of 1 μm² s^–1 as we have estimated for a diffusible signal, Unpaired, in the same biological system ^[29], and near the ranges determined experimentally for diffusible morphogens in Drosophila ^{[31,32,33,34]}. Different values of k were also tested to find their impact on the time to the 1 pM mark at the anterior. The parameter k is the degradation of chemoattractant in the extracellular space. The values tested were 1 × 10^–2 s^–1 to 1 × 10^–6 s^–1 on a logarithmic scale. With k values of 1 × 10^–2 s^–1 and 1 × 10^–3 s^–1, the time for the chemoattractant to diffuse across the egg chamber surpassed the five hour mark. This suggests that k values of 1 × 10^–2 s^–1 and 1 × 10^–3 s^–1 may be too high. For k sufficiently small we see little change in time to reach the 1 pM threshold. A value of k near 1 × 10^–4 s^–1 seems to be a reasonable order of magnitude for the degradation rate of an extracellular signaling molecule ^{[29,32,33,34]}. The next parameter we tested is ϕ. Values from 1 × 10^–1 μm s^–1 to 1 × 10^–5 μm s^–1 were tested for the uptake of chemoattractant into the nurse cells, ϕ on a logarithmic scale. We found that for a value of ϕ = 0.1 μm s^–1, the chemoattractant did not reach 1 pM at the anterior of the egg chamber in under five hours, suggesting that ϕ should be less than 0.1 μm s^–1. Similar to our results for k, small values of f caused minimal change in the time to reach threshold. The secretion constant σ is the amount of secreted chemoattractant per unit of space in a given time. The values of σ that we tested ranged from 1 pM μm s^–1 to 10,000 pM μm s^–1 on a logarithmic scale. We found that the higher values of σ resulted in less time for the point in the anterior of the egg chamber to reach a 1 pM concentration, as would be expected. Changing this variable resulted in a large range of times for the chemoattractant to reach the end, as can be seen in Table 1. However, it is important to note that σ did not make the chemoattractant diffuse faster but rather amplified the concentration of the chemoattractant participating so that it rose above the 1 pM anterior level in fewer seconds; σ does not change the diffusion but rather it changes the magnitude of the chemoattractant influx.

Figure 4 shows an internal view of the distribution of chemoattractant at the time point at which the signal reached threshold at the anterior end of the egg chamber. The nurse cell interiors are not accessible to diffusing chemoattractant. In this case, the chemoattractant is being secreted from the full face of the oocyte, as well as the epithelium where x > 300 μm. The model parameters are D = 1 μm² s^–1, k = 1 × 10^–4 s^–1, f = 0 μm s^–1, and σ = 100 pM μm s^–1. It is noteworthy that, while the migratory path is anterior to posterior (left to right) through the center, the areas between nurse cells accumulate significant levels of chemoattractant in the dorsal-ventral axis (up and down), for example at x = 200 and x = 300 (Figure 3, bottom panel), which may impact directional migration. Some studies in vivo have shown the migratory cells alter their behaviors at these points, for example switching positions with each other. However, the gradient is steeper toward the posterior, which may explain how the cells navigate correctly at these locations. To test the impact of chemoattractant distribution along the dorsal-ventral axis and determine if this causes the rotational movements, it would be interesting to mis-express chemoattractant to the ventral side and observe changes in migratory cell behaviors in vivo.

---

3.2. Equilibrated chemoattractant distribution and gradient occurs within two to three hours

When the concentration of chemoattractant at each point reaches a level where it will consistently remain, it has achieved a steady state. Similarly, when the gradient of the concentration no longer changes, it has reached a steady state. To estimate our rate of change, we define an error estimate, find a vector of the concentration or gradient at every element at a particular time step, subtract the same vector for the previous time step, and take the Euclidean norm. When this value goes below a certain cutoff, we decide that the chemoattractant has reached a steady state. We used a cutoff of 0.0008 (Table 2).

To calculate the direction of the gradient, we first normalized it. Then when the average difference of the normalized gradient between timesteps fell below a cutoff (e.g., 10^–2), we recorded that time. As with the concentration, the measure of changing gradient direction decreases over time, meaning the direction of the gradient becomes more stable (Table 3). With our definition of the steady state of the concentration and of the directional gradient, the directional gradient comes to a steady state typically more quickly than the concentration does; the direction often can be established more quickly.

---

3.3. Varying parameters impacts the average amount of chemoattractant in the egg chamber extracellular space

Multiple cell types in the egg chamber, including the oocyte and follicle cells, express the genes encoding EGF-like chemoattractants, but it is not clear exactly which cells process and secrete the proteins ^[24,25,40]. We considered that the epithelial cells along with the surface of the oocyte may be secreting the chemoattractant, so we tested how changing the region of epithelial cells secreting and the total amount of chemoattractant secreted would affect the level of concentration in the egg chamber after 5 hours. We were searching for a combination that would produce an average chemoattractant concentration of approximately 4500 pM, an in vitro experimental activating value for other secreted signaling molecules ^[35], and in the range of chemotactic-activating concentrations for EGF ligands ^{[36,37,38,39]}. With this goal, we changed the approximate number of epithelial cells that were secreting the chemoattractant by altering which locations of the epithelium were actively secreting during certain trials. Since the epithelial cells condense toward the oocyte while the gradient is being generated, secretion from more dense, posterior regions of the epithelium may effectively be higher. Considering that the coordinates of the egg chamber in the x direction go from 0 μm to 409.5 μm, we used epithelium locations including the oocyte and starting at 200 μm, 250 μm, 300 μm, 350 μm, and just the oocyte (at 409.5 μm). Then, for each of these chosen active areas, we changed the value of σ, the flux of chemoattractant entering the egg chamber, to determine the average chemoattractant concentration in the egg chamber after 5 hours under each condition (Table 4).

Many conclusions could be drawn from such a test. As expected, keeping the area of secreting epithelium cells the same, as we increased σ, the concentration in the egg chamber increased proportionally. Keeping σ the same and increasing the amount of the epithelium secreting also caused the final average chemoattractant concentration in the extracellular space of the egg chamber to increase. We then tested the same σ values and locations of secretion with two different values of k and φ, the extracellular degradation and nurse cell uptake of chemoattractant respectively. When we increased k and φ by a factor of 10, changing them from 1 × 10^–4 s^–1 to 1 × 10^–3 s^–1 and from 1 × 10^–4 μm s^–1 to 1 × 10^–3 μm s^–1 respectively, we observed the concentration average decrease for corresponding combinations of secretion location and σ value. This is logical because more of the chemoattractant is taken up by the nurse cells. However, it is interesting to note that it cuts the average concentration almost exactly in half.

Looking at Figure 5, each pair of adjacent bars, starting with any tall one and looking at the shorter one directly to its right, shows that increasing k and φ by a factor of 10 and doubling σ, but maintaining all other parameters will cause the average concentration to be cut in half. Therefore, when thinking about the value of σ needed to produce an average concentration level of approximately 4500 pM after 5 hours of diffusion, maintaining the same secretion location and increasing k by a factor of 10 will necessitate doubling the value of σ in order to reach this realistic chemoattractant value. So, for a measured average concentration of chemoattractant, these tables can be used to find appropriate secretion values. More experimental tests will be needed to determine where the chemoattractants originate, which can be used in combination with this result to estimate values.

---

3.4. Large chemoattractant difference across egg chamber precludes a linear migratory response to gradient

A gradient of chemoattractant is understood to be necessary for cell migration. The mechanism by which the cells interpret directional cues from the gradient is unclear, although several models have been proposed (for example ^[4,17,30,41]). We considered the possibility that the cluster responds linearly to the magnitude of the gradient. To determine if this is plausible, we tested the magnitude of the x-component of the gradient along with the chemoattractant concentration at points near the anterior and posterior ends of the egg chamber. We chose a point near the posterior end (Back: x = 400, y = 150, z = 150) and a point near the anterior end (Front: x = 20, y = 150, z = 150) (as in Figure 2) and found the ratio of the gradient and concentration (Tables 5 and 6, respectively).

When k and ϕ had higher values, such as 0.01 μm s^–1 and 0.01 s^–1 or 0.1 μm s^–1 and 0.1 s^–1, there was a greater difference in concentration between the anterior and posterior because more chemoattractant would deteriorate before reaching the end. When k = 1 × 10^–1 s^–1, there was too little concentration at the anterior end to give a reliable value. In some nearby locations the y or z component of the gradient pointed in opposing directions (Figure 6), which could cause disruption during directional migration due to competing migratory forces experienced on different sides of the migrating cluster of cells. The concentration decreased everywhere, but less so near the oocyte. A very high D, such as 1000 μm² s^–1, gave a concentration that was lower near the oocyte, but higher in the anterior than with lower D because the chemoattractant was able to diffuse more freely. Changing s causes no difference in the ratio of either concentration or gradient because it is only a scaling factor. Using a larger area of secretion makes the ratio decrease. In general, changing concentration ratios corresponded with changing gradient ratios.

No matter how we changed the parameters within the biologically plausible ranges, every simulation resulted in gradients at opposite ends of the egg chamber whose value differed by at least one order of magnitude. We concluded that the migratory force cannot scale linearly to the magnitude of the gradient because in vivo experiments show a less than a two-fold difference in migration speed ^[21,22,42]. This distribution agrees with steady state, one-dimensional estimates of a more exponential chemoattractant distribution described in ^[27]. It may be that there is a low saturation point where the force stops increasing linearly, the force is very nonlinearly related, or that the migratory force only depends on the direction of the gradient and not the magnitude.

Next, we analyzed the chemoattractant concentration and the x, y and z-components of the gradient over time, along the center of the egg chamber, which is the migratory path for border cells (Figure 7). We observed that a gradient of chemoattractant is established quickly, within the first 30 minutes of the 5-hour timespan. Interestingly, separation of the spatial components of the gradient revealed major declines in the x-component near nurse cell junctions and near the oocyte, with reversals in the y and z-components nearby, particularly near x = 300 μm and x = 390 μm. These results suggest a significant impact of interstitial spaces between nurse cells on the shape of the chemoattractant gradient. Some studies in vivo have shown the migratory cells alter their behaviors at these points, for example switching positions with each other (or tumbling) ^{[21,22,26,27]}. The weak reversals of gradient could help to explain this unexpected behavior.

---

3.5. Egg chamber without nurse cell packing exhibits a less steep gradient

An important question about cell migration is how the nurse cells interact with the diffusion of the chemoattractant. In our model, the nurse cells either act as an impenetrable barrier and completely impede the chemoattractant, when ϕ = 0, or uptake some level of chemoattractant when ϕ is nonzero. Our model leaves space (approximately 19.5 μm, the diameter of 2.8 IMCs and the approximate diameter of the migratory border cell cluster) between nurse cells for the chemoattractant to travel, but in reality, it can be inferred that the interstitial pathways are likely to be narrower. We also considered what would happen to the established gradient if the chemoattractant could travel freely through the nurse cells or be absorbed by them as a part of a sink in the egg chamber. To get an idea of the effect of nurse cells blocking chemoattractant, we tested the concentration and x-component of the chemoattractant at different points (at x = 20 μm; 100 μm; 200 μm; 300 μm, 400 μm) along the central corridor after a simulation was completed on one geometry with nurse cells and one without them. We found that the chemoattractant is more evenly dispersed without border and nurse cells (which potentially mimics the effect of chemoattractant moving freely through entirely permeable nurse cells) (Table 7). There is a greater disparity between the concentration values at the posterior and anterior ends of the egg chamber after the diffusion with the nurse cells impeding chemoattractant; the same was found for the x component of the gradient. We can also observe that the nurse cells in the geometry cause the gradient to be stronger than what it would be had the nurse cells not been present (Table 7). A similar effect on the gradient was observed for models in which the egg chamber domain was half of the size (200 μm, Supplementary Table 1 and 2). We believe this result arises not because of dilution effects but because the nurse cells provide a barrier that prevents the chemoattractant from freely spreading farther away from the oocyte face. The results of this experiment imply that the narrower the domain in which the chemoattractant can disperse, the larger the range of the gradient and concentration in the egg chamber.

---

3.6. Cell migration kinetics are altered by chemoattractant distribution model

To obtain a sense of how of migratory cells respond to the chemoattractant distribution using our cell-packed geometry, we used a previously-described collective cell migration model ^[26] and input resultant chemoattractant concentration gradients derived from the chemoattractant model. Figure 8 (left) shows the starting position of the group of 6 migratory cells in the anterior domain of an egg chamber in three dimensions, and (right) shows the large nurse cells differentiated by colored IMCs. Figure 9 shows the early phase of migration with the chemoattractant in a cell-packed model or a unit gradient. When signal distribution is modeled in a realistic tissue geometry, the cells migrate more slowly. In the early phase, the gradient is shallow and cells move slowly; later, the cells' speed increases exponentially, reflecting the exponential nature of the gradient.

---

4. Conclusions

Many biological systems rely on diffusible, extracellular signals to convey diverse kinds of information to cells. For example, migratory cell types respond to gradients of chemoattractants and chemorepellants. To understand this response, it is critical to determine the spatial distribution of signals within a complex tissue or biological system. However, these signals can be present at levels that are not easily detected, and methods to measure them spatio-temporally in vivo are inadequate. Mathematical modeling provides a way to gain insight into the possible distribution of these molecules within the complex geometry of animal tissues.

We developed a mathematical model of the well-described Drosophila egg chamber to examine possible chemoattractant distribution in this context. A set of motile cells called border cells chemotax towards several cues produced from and around the oocyte, but the exact distribution of the cues is unclear. Large nurse cells line the path of the border cells and likely impede movement of the chemoattractant signals (Figure 1). Using estimates for chemoattractant diffusion rate, decay, secretion rate, and uptake by other cells, we determined reasonable parameters that would allow the chemoattractant to reach the border cells in the appropriate timeframe to guide their movement. Previous work postulated that chemoattractants move uniformly through all cells of the egg chamber to generate a linear or exponential gradient ^[26,27]. We contrasted that model with one in which signals did not move through the nurse cells, but instead were blocked by them.

Notably, the physical constraint of the nurse cells dramatically increases the local concentration of chemoattractants along the migratory route and decreases the time at which the gradient sets up across the tissue, compared to a case in which nurse cell packing was not considered a hindrance. In the model, due to computational limitations with meshing inflated nurse cells, we overestimated the space between nurse cells and nurse cells and epithelium, however, a small deflation of the nurse cells only affected the time to reach the 1 pM threshold at the anterior end of the cell for sensitive ϕ_eff = 0.36 taking more time, about 4.95 hours compared to 3.6 hours. The time to threshold while changing both D_eff and σ_eff were unaffected by this slight deflation. In general, the smaller gaps between the cells seen in vivo would result in even more dramatic differences in chemoattractant distribution (compared to non-hindrance, see Supplementary Tables 1 and 2). This idea has important biological implications since cellular obstacles that alter extracellular signal distribution would effectively concentrate the signals, allowing development to proceed more quickly, and would require less signal to be produced, conserving cellular energy.

The idea that tissue geometry impacts signal concentration is made more important by recent studies that suggest ligand concentrations have key impacts on signal transduction. In both mammals and Drosophila, multiple ligands are capable of binding to the guidance receptor tyrosine kinases. Recent studies suggest that different ligands at the same concentration produce more similar downstream effects than the same ligand at different concentrations ^[36,43] similar to how morphogens are thought to act in a concentration dependent way ^[2,3,32]. Thus, it is essential to consider how tissue shape may alter the perceived concentrations of signaling molecules.

Interestingly, live imaging studies show that migratory border cells alter their movements at the points where multiple nurse cell meet, exhibiting "tumbling" behaviors ^[21,22]. This could in part be explained by physical constraints that change available traction points for the migratory cell, as indicated by Stonko et al. ^[26], and other possible molecular explanations have been proposed ^[21,22]. The model presented here suggests that changes in the shape/steepness of the gradient could also influence this behavior, as we detect drops in the gradient at the junctions between nurse cell boundaries. Most cell migration models presume the cells are responding to a linear or exponential, steadily-increasing chemoattractant gradients ^{[4,6,7,41,44]}, but in vivo it is likely that non-uniform landscapes of cellular architecture or ECM could produce differences in signal distribution that need to be considered more carefully. This could be important to determine for animal development, and also immune response and tumor metastasis ^[5,45].

In pilot tests, collective cell migration simulations show that chemoattractant distribution impacted by cell boundaries also impacts the kinetics of cell movement. Figure 9 shows the early phase of migration with the chemoattractant in a cell-packed model or a unit gradient. When signal distribution is modeled in the realistic tissue geometry, the cells migrate more slowly overall. The exponential nature of the gradient is realized in the migration with speed increasing as the cluster nears the oocyte. For time-dependent developmental or immunological migrations, this difference in response could have a dramatic phenotypic output. Additional in depth simulations and analysis will be needed to determine the details of individual cellular responses to the distribution, and to ascertain if the slower movement is due to the more shallow gradient or variations in gradient magnitude and direction.

Currently, it is technically challenging to detect diffusible molecules accurately in vivo. Prior work has shown that diffusible signaling molecules are clearly restricted in tissues, moving much more slowly than in aqueous solution. Experimentally, the measured values vary widely, from 0.1 to 100 μm²/s, in part depending on the methods used ^{[1,33,34,46,47]}. Chemoattractants have not been examined in vivo in living egg chambers, and for EGF ligands, the predominant sources of the secretion is unclear, which is why we modeled different possibilities for domains of release. With either secretion from the oocyte only or from the oocyte and follicle cells, while keeping the other parameters the same, we found reasonable timeframes to activate cell migration. The impact of the different origins of signal on migration will be interesting to explore further.

The egg chamber provides a unique context in which to study the impact of cell packing on signal transduction. This work is facilitated by the diversity of cells in this tissue, and the ease with which Drosophila can be genetically manipulated and imaged live. While this tissue may be somewhat unique in structure, we believe interstitial gaps may be prevalent in animal tissues but difficult to discern with common imaging techniques. Thus, the impact of the cytoarchitecture on signaling should be more widely considered and modeled.

---

Acknowledgements

These results were obtained with support from REU Site: Interdisciplinary Program in High Performance Computing (hpcreu.umbc.edu) in the Department of Mathematics and Statistics at the University of Maryland, Baltimore County (UMBC). This program is funded by the National Science Foundation (NSF-1460652), the National Security Agency (NSA), and the Department of Defense (DOD), with additional support from UMBC, the Department of Mathematics and Statistics, the Center for Interdisciplinary Research and Consulting (CIRC), and the UMBC High Performance Computing Facility (HPCF). HPCF is supported by the U.S. National Science Foundation through the MRI program (grant nos. CNS-0821258 and CNS-1228778) and the SCREMS program (grant no. DMS-0821311), with additional substantial support from UMBC. J.C. was supported, in part, by the UMBC National Security Agency (NSA) Scholars Program through a contract with the NSA. M.S. was supported by UMBC. M.S.G. and B.E.P. were supported in part by grant NSF-IOS-1656550. We thank Dr. Daniel Lobo for helpful comments on the manuscript.

---

Conflict of interest

The authors declare no conflict of interest in this paper.

---