A New Chemotactic Mechanism Governs Long-Range Angiogenesis Induced by Patching an Arterial Graft into a Vein

Chemotaxis, the migration of cells in response to chemical stimulus, is an important concept in the angiogenesis model. In most angiogenesis models, chemotaxis is defined as the migration of a sprout tip in response to the upgradient of the VEGF (vascular endothelial growth factor). However, we found that angiogenesis induced by performing arterial patch grafting on rabbits occurred under the decreasing VEGFA gradient. Data show that the VEGFA concentration peaked at approximately 0.3 to 0.5 cm away from the arterial patch and decreased as the measurement approaches the patch. We also observed that the new blood vessels formed are twisted and congested in some areas, in a distinguishable manner from non-pathological blood vessels. To explain these observations, we developed a mathematical model and compared the results from numerical simulations with the experimental data. We introduced a new chemotactic velocity using the temporal change in the chemoattractant gradient to govern the sprout tip migration. We performed a hybrid simulation to illustrate the growth of new vessels. Results indicated the speed of growth of new vessels oscillated before reaching the periphery of the arterial patch. Crowded and congested blood vessel formation was observed during numerical simulations. Thus, our numerical simulation results agreed with the experimental data.


Introduction
Angiogenesis is the formation of new blood vessels that occurs under the control of chemical signals in the body. It is responsible for a wide variety of physio and pathological processes and the progression of diseases such as cancer, cardiomyocytes, or arteriovenous malformations (AVMs). Vascular endothelial growth factor (VEGF) is one of the chemical signals that bind to receptors on the surface of endothelial cells [1]. When VEGF binds to the surface receptors, endothelial cells become activated and detached from the blood vessel and subsequently migrate to create new blood vessels. The migration of new blood vessels is regulated by the spatial gradient of the signal. This process is called chemotaxis and the signal is called a chemoattractant. The endothelial cell at the front of the sprout tip senses the VEGF gradient to lead the direction of migration. This migration is then followed by other endothelial cells that then create a path line of a new blood vessel. Typically, in angiogenesis, chemical signals are released at the time and place where they are required. Therefore, new blood vessel formation is well-balanced in processes such as growth and healing. However, in the state referred to as the angiogenic switch [2], the chemical signals can occasionally become unbalanced and irregular, causing abnormalities in the new blood vessel formation, leading to vascular anomalies. For example, arteriovenous fistulas (AVFs) and arteriovenous malformations (AVMs).
Arteriovenous fistulas (AVFs) are the anomalies in which arteries and veins are directly connected, bypassing the capillary network that results in the shunting of blood between the two. AVFs are distinct from AVMs and are congenital lesions that arise due to failure of embryonic or fetal vascular differentiation and are also related to noniatrogenic trauma [3,4]. Usually, they increase in size if untreated. As for arteriovenous malformation (AVM), it is an abnormal angiogenesis consisting of feeder artery and drainage vein vessels which is abnormal, complex, and called "Nidus". AVMs probably have high flow compared to AVFs, higher resistance, and may remain asymptotic [3]. AVMs can occur in organs of the human body such as the brain, lungs, intestinal tract, and the spinal dura mater [5][6][7][8]. So far, the cause of AVMs is not entirely clear. In particular, the process of AVF formation, which is a crucial step in AVM, pathogenesis is not yet understood. Moreover, the role of angiogenic factors, including the VEGF, in AVMs is still unknown owing to the lack of spontaneous in vivo AVM models. In fact, the lack of in vivo AVM models affected the development of therapeutic medications for AVMs [9].
Recently, ref. [10] have successfully induced AVF by patching an arterial graft into a rabbit's vein. The new blood vessels were distinguishable from the pre-existing linear blood vessels owing to strong flexion and significant meandering. In this study, we extend the study in [10] by analyzing the VEGF concentration during AVF formation in an arterial patch procedure. The VEGF concentration was measured at a certain period of time. We observed that the VEGF concentration oscillated and gradually attenuated with time. Based on this data, we developed an angiogenesis model by introducing a new chemotactic velocity that governs new blood vessel formation under certain VEGF concentration changes.
In this study, we took a hybrid approach and modified Anderson and Chaplain's model [23]. In their model, it was assumed that the migration of sprout tip is influenced by random motility, chemotaxis in response to signal gradients, and haptotaxis in response to fibronectin gradients. As shown in the following equation, three factors contribute to the endothelial cell flux J n : They described the chemotactic flux J chemo as where n, c, and χ denote the sprout tip density, chemoattractant concentration, and chemotactic coefficient, respectively. Generally, the conservation equation for the sprout tip density is in the form of In Section 2.5, we expressed spatiotemporal molecular concentrations as partial differential equations and individual cell movements as biased random walks. We found that our simulations were consistent with the experimental observations.

VEGF Distribution Model
VEGF is the key mediator for angiogenesis and has been discovered to be a permeabilityenhancing agent, leading to the disruption of intercellular contacts and increase of permeability [24]. The VEGF family of genes consists of at least seven various subtypes such as VEGFA, VEGFB and VEGFC, VEGFD, PIGF, VEGFE (Orf-VEGF), and trimeresurus flavoviridis svVEGF. However, VEGFA and its receptors, vascular endothelial growth factor receptor (VEGFR), including VEGFR-1 and VEGFR-2, play a major role in pathological angiogenesis. VEGFA has a variety of functions such as vascular permeability activity and stimulation of cell migration. Here, VEGF is also known as VEGFA [25].

VEGFA Quantification in In Vivo Angiogenesis Model
The in vivo angiogenesis-inducing model in [10] sutures an arterial graft into the wall of the left common jugular vein in male rabbits for the initiation of angiogenesis. Using this model, we quantified the spatiotemporal VEGFA concentration with an antibodybased measurement method called ELISA. Detailed experimental protocols are given in Section 4.2. In this paper, we consider the VEGFA concentration as the representative of the VEGF concentration and use it as the basis for our mathematical modeling.
The VEGFA concentration at the arterial patch graft was low and the level was equal to that of the serum. The highest level of VEGFA concentration was observed near the pre-existing artery (F4/5). At any time point, VEGFA concentration was higher in the pre-existing artery (F4/5) than in the arterial patch graft (Fp) (see Figure 1b). VEGFA concentration also revealed a repeated increase and decrease in addition to gradual attenuation, where the approximate curve shows a sine wave that attenuates gradually (see Equation (4) below). In this angiogenesis-inducing model, angiogenesis originates from the second branch of the left subclavian artery, heading for the arterial graft in the common jugular vein. From these observations, it is presumed that VEGFA is primarily released from the area where new blood vessels sprout and grow. In other words, blood vessels migrate toward the VEGFA with low gradient concentration.
(a)  , and F7) that have been taken to measure the VEGFA concentration. The arterial patch graft is placed at Fp and the pre-existing artery is assumed to exist at around F4/F5 areas. (b) VEGF concentration from experimental data showing that the VEGF concentration was higher in the pre-existing area (F4/5) than in the patch area (Fp). (c) The measured VEGF concentrations (pg/mL) at the pre-existing vessel change with time. This indicates that there were three rabbit cases for each time point. We plotted each case on the graph, and its fitted graph was derived. The curve was obtained by numerically fitting the data to Equation (4): a = 46.2525, b = 0.2131, σ = 0.1548, γ = 0.3485.

VEGFA Distribution Equation Based on In Vivo Quantification Data
Let c(x, t) be the distribution of VEGFA concentration at position x ∈ [0, 1] or [0, 1] 2 and time 0 ≤ t ≤ T. Here, the VEGFA distribution equation is based on the VEGFA quantification of the arterial patch grafting in Section 2.2. Based on the quantification above, we assume the following conditions for the VEGFA distribution:

1.
The VEGFA concentration near the pre-existing vessel is always higher than that near the arterial patch ( Figure 1a).

2.
The VEGFA concentration changes according to a sine wave that attenuates gradually, as in Figure 1b.
Define the distribution of VEGFA concentration c(x, t) as follows: where a, b, σ, γ, and are constants. Here, r p is the distance function from the outer side of the arterial patch, which is defined as where (x p , y p ) is the location of the patch center. We assumed that the arterial patch has a radius of 0.1 (dimensionless) and was placed on the left side of the boundary (i.e., x p = 0). Figure 2b shows the modeled time course of VEGFA concentration c(x, t) near the pre-existing artery (x = 1). Figure 2c shows the initial VEGFA distribution c(x, t) (t = 0) at domain (x, y) ∈ [0, 1] 2 . We used this VEGFA profile for a numerical simulation of the angiogenesis model.    (4). The VEGFA concentration is always higher near the pre-existing blood vessels other than in the other area. This illustration is based on the experimental data in Section 2.2. (b) VEGFA concentration profile using Equations (4) and (5) at F4/F5 by setting parameters as follows: x, y ∈ (0, 1), a = 1, b = 0.39, σ = 0.23/2, γ = 0.15, = 2.9, x p = 0, y p = 0.5, x = 1, and y = 0.5. This VEGFA concentration changes graph has a qualitatively similar profile as that in vivo. The VEGFA concentration changes graph has a similar profile as in vivo experiments (see Figure 1b). (c) VEGF distribution of Equations (4) and (5) on domain x, y ∈ (0, 1) at t = 0.

New Blood Vessel Length Measurement in In Vivo Angiogenesis Model
Here, we measured the new blood vessels with anomalous vessel walls sprouting from pre-existing arterioles (see Figure 3a). As shown in Figure 3b, the lengths of the new blood vessels increased gradually until day 3, while only a small change of length was observed from day 3 to 5. On day 6, the vessels grew again, and the growth continued until day 9. The new blood vessels did not grow much further after day 9. In this model, the arteriovenous fistula was formed between days 10 to 14 [10]. Therefore, we speculate that the angiogenesis process stopped once the new blood vessels that sprouted from the arterioles reached the arterial graft on day 9. (a)

Sprout Tip Distribution Equation with New Chemotactic Velocity
We assumed that the migration of the sprout tip was influenced by random motility and chemotaxis in response to temporal changes in the VEGFA gradient, rather than the gradient itself. To derive the sprout tip distribution governed by these two factors, we considered the following mass conservation of sprout tip density n = n(x, t): where d denotes the random motility constant and v m = v m (x, t) is the transport velocity [23]. Equation (6) is divided into two terms. The first term explains random motility. The second term explains the migration of the sprout tip under the influence of a flow with velocity v m . Partly inspired by the analysis of the LEGI (local excitation, global inhibition) where Here, β is constant and c t in (8) is the derivative of our VEGFA distributions (4) and (5). The application of (8) was also inspired by the comparison of the growth of new blood vessels ( Figure 3b) and the approximation curve of VEGFA concentrations (Figure 1c). The growth speed of new blood vessels was relatively fast from day 0 to 2 and from day 5 to 8, which seems to coincide with the timing when the VEGF concentrations increased near the pre-existing vessel. However, we mentioned that we cannot decisively say the above two must be correlated, owing to relatively high variability in the vessel growth data (Figure 3b).

Discussion
Here we used the temporal change in the chemoattractant gradient as the new chemotactic mechanism to model the arterial patch-induced angiogenesis. The temporal changes in VEGF concentration affected the angiogenesis process. This was probably caused by the hypoxia state around the arterial patch, which affected the VEGF attenuates status to stimulate the sprout tip of the endothelial cells. Our angiogenesis model successfully predicted new blood vessel growth that fitted the experimental data. During the numerical simulation, we successfully obtained twisted and congested formations that were formed from t = 2 to t = 7. This formation was pathological blood vessels characteristic.
We compared the length of the new blood vessel at each time during the numerical simulation, from t = 1 to t = 14 with the experimental results [10], as shown in Figure 5. Curves in different colors show different simulation results under the same parameter sets (see Figure 6b). Although we used dimensionless parameters for the numerical simulation, the simulation results showed patterns that were similar to those from in vivo measurements (see Figure 6). Because of the lack of experimental data and little information, we can see a differentiation value in error for the mean of each new blood vessel's growth from in vivo. The growth of the new blood vessels was initially fast before decelerating at some time period. Subsequently, the new blood vessels regrow until reaching the arterial patch. Additionally, we also ran our simulations with different parameters since our method also depends on spatial grid step discretization (see Figure S1).

Derivation of Numerical Scheme
To understand our angiogenesis model, we solve systems (4)-(8) numerically and perform a simulation of the new blood vessel growth using a hybrid technique. This technique was invented by Anderson and Chaplain to simulate a tumor-induced angiogenesis model [23]. We herein generalize the hybrid technique such that the scheme can be used to simulate our systems (4)-(8) for any transport velocity v m .
We consider the following system: and define the domain of system (9) as Ω = (0, 1) × (0, 1). The following numerical scheme is known to ensure both positivity and mass conservation of the solution [26]. The domain Ω is partitioned into two types of lattices, i.e., main and sub: where h = 1/N (see Figure 4a). For the time step, we renew the time increment, τ k , at every step k.
where τ k will be defined later. i,j+1 , and P k,4 i,j−1 to P k,1 i,j , P k,2 i,j , P k,3 i,j , and P k,4 i,j , so that each of these P on (i, j) can be assumed as the probability of a sprout tip to migrate to the left, right, up, and bottom, respectively. Additionally, P k,0 i,j is regarded as the probability of a sprout tip to stay on (i, j).
In these lattices, we assume that n and v m are approximated on the main lattice. Additionally, we assume that f , G, and F are approximated on the sub-lattice: Here, (v m ) k i,j = (v m ) k i,j,1 , (v m ) k i,j,2 is calculated on the x and y directions using the average of four points surrounding every (1 ≤ i, j ≤ N).
where f k i,j = f (x k i ), and f is calculated on the sub-lattices. For f (x, t) in (8), (c t ) k i,j can be calculated by taking the time derivative of (8) on the sub-lattices, for all (1 ≤ i, j ≤ N) as the direction setting of the upwind side. Next, we consider the sub-lattice where G k i, is defined. In the x direction, we assume that n k i,j and n k i+1,j are brought to G k i,j by flows (v m ) k,+ i,j,1 and (v m ) k,− i+1,j,1 , respectively. We may also assume that n k i,j and n k i,j+1 are brought to G k i,j by flows (v m ) k,+ i,j,2 and (v m ) k,− i,j+1,2 , respectively, in the y direction; see Figure 4b. Thus, we write G k i,j = (G k i,j,1 , G k i,j,2 ) as follows: Therefore, F k i,j = (F k i,j,1 , F k i,j,2 ) can be approximated by the following scheme.
By the zero-flux in Equation (5) of a system (9), we have the following conditions: Hence, n k i,j is approximated by the following numerical scheme: We substituted (13), (15), and (16) into (18), which results in the following scheme: n k+1 i,j = n k i,j P k,0 i,j + n k i+1,j P k+1 i+1,j + n k i−1,j P k,2 i−1,j + n k i,j+1 P k,3 i,j+1 + n k i,j−1 P k,4 i,j−1 , where for (2 ≤ i, j, ≤ N − 1). The same procedure also applied for i, j = 1 . . . N to obtain the scheme at the boundary. From the scheme above, the values of P k,0 i,j , P k,1 i+1,j , P k,2 i−1 , P k,3 i,j+1 , and P k,4 i,j−1 are the weights from each point of (i, j)'s neighbors to approximate the value of n k+1 i,j . Assume five P values as the probability of a cell remaining at (i, j), migrating to the left (i − 1, j), to the right (i + 1, j), to the downside (i, j − 1), and to the upside (i, j + 1), respectively, if we shift P k,1 i+1,j , P k,2 i−1,j , P k,3 i,j+1 , and P k,4 i,j−1 to P k,1 i,j , P k,2 i,j , P k,3 i,j , and P k,4 i,j , respectively, see Figure 4d.
During angiogenesis, a new sprout tip can be generated from the existing sprout tips. Additionally, the sprout tips can form loops by connecting with other sprout tips or the existing new blood vessels. These two events are called branching and anastomosis, respectively [27]. We adopted the branching and anastomosis rules from the discrete mathematical model of tumor-induced angiogenesis [23]. Namely, for the branching to occur, the lifetime of sprouts must exceed the threshold age T branching , and the branching probability is given by with a modification that the branching probability is applied only when the temporal change in the VEGFA concentration is positive, which follows (8). When anastomosis occurs, only one sprout is allowed to continue to grow. As detailed in [26], this numerical scheme is stable when the time increment τ k is updated by the following criterion: where and P k,l i,j ≥ 0 (l = 0, 1, 2, 3, 4), which confirms the probability assumptions of Ps from schemes (19)-(24) is reasonable.

Numerical Simulation Results
In this section, the domain setting and parameter values used in the simulations are depicted. All the parameters are dimensionless. As discussed in Section 2.3, we assumed that an arterial patch of radius 0.1 (dimensionless) was placed in the center of the left boundary. Thus, we set r p = x 2 + (y − 0.5) 2 . Initially, we placed 20 sprout tips randomly at the right boundary. This position can be considered to be that of the pre-existing vessel (see Figures 1a and 2a) by setting the length of domain-x as the actual in vivo length.
We show our numerical simulation results of the angiogenesis model (4)-(8) in Figures 5, 6 and 7, respectively. Figure 8 shows the initial state of the arterial patch on the left side of the boundary and pre-existing vessel on the right side of the boundary at t = 0. The VEGF distribution is higher near the pre-existing vessel and lower towards the arterial patch. We simulated for T = 0 until T = 14 and tracked down the growth of the new blood vessels at t = 1, 2, 7, 10, 11, and 14. In these figures, the yellow to orange color in the background of each picture represents the VEGFA concentration (see also Equation (4) and Figure 8). In Figure 5, we see that at t = 1 and t = 2 the VEGFA concentration is higher (orange), but it gets lower (yellow) at later times. Since the VEGF concentration attenuates gradually (see Figure 9a) over time followed by the sprout tips of the newly generated blood vessels, we can see changes in color from orange to yellow and vice versa.
We looked at the blood vessel growth in Figure 5. At t = 2, the growth of new blood vessels started to decelerate and the sprout tips moved randomly. The random movement of the sprout tips heightened the probability of the blood vessels meeting other tips which caused the congested and twisted formation. This creates an AVM formation such as the one in Figure 5c. Thus, many sprout tips died owing to the anastomosis. The fewer the sprout tips generated initially, the lower the probability of anastomosis. In fact, there were two sprouts at most (from 20 sprout tips generated initially) that survived and reached the arterial patch resulting in vascularization (see Figure 5d) during the numerical simulation. Figure 6 shows temporal changes in the VEGF gradient causing the graphs to have a high and almost linear slope between in vivo experiments and our simulations. We ran a few simulations and calculated the vessel length in millimeters (mm). We can see that the two graphs are almost similar. Both graphs show that at day 2, the growth of the sprout tip started to decline slowly but regrow faster on day 7 until it reached the periphery of the arterial patch on day 10.
Data for each new blood vessels length from simulations and in vivo are presented in the graph (see Figure 7) of mean ± standard deviation (SD), 1.04571 ± 0.47451 and 0.80 ± 0.22234, respectively. Statistical analysis was evaluated by Student's t-test and the p value was < 0.05. Figure 7b shows that the result was statistically significant for each day. As shown in Figure 7b, we can see from day 2 until day 7 the growth speed of blood vessels becomes slow and almost constant and the speed becomes faster and higher after day 7 and become constant again starting from day 10. The attenuated status of VEGF also affected the growth speed of the blood vessels. The p value in vivo was >0.05, which is statistically not significant (see Figure 7a). However, our simulation results show similar patterns to those in vivo experiments ( Figure 5).   Figure 9. Comparison of the new blood vessel growth and the VEGFA gradient changes. When the VEGFA concentration is decreasing over time, the growth of the new vessels become slower. We consider this growth of the new blood vessels at the decreasing rate of VEGFA to be called the rest period, R p , and at the increasing rate of VEGFA is called the growth period, G p . (a) VEGFA concentration changing over time (days) and attenuates gradually. (b) Measurement of the new blood vessel in days.

Quantification of Vegfa in In Vivo Angiogenesis Model
Three rabbits at each time point (1, 3, 6, 12, 18 h, 1-7 days; different rabbits each time) after the grafting was euthanized, and the arterial graft and fatty tissue were collected and then divided into approximately eight areas, as shown schematically in Figure 1a. The total protein was extracted from the tissues by the Minute Detergent-Free Total Protein Extraction Kit obtained from Invent Biotechnologies, Inc. (Minneapolis, MN, USA). The level of VEGFA concentration in the tissue was measured by Rabbit VEGFA ELISA kits purchased from Cusabio Biotech Co., Ltd. (Wuhan, China). The concentration of VEGFA in each sample was calculated based on a recombinant protein-based standard curve (Supplementary Materials Figure S4). The Multiskan JX spectrophotometer (Thermo Fisher Scientific Inc., Yokohama, Japan) was used to measure the absorbance. All assays were performed in duplicate, and the mean values of the results were used in the analyses.

Measurement of New Blood Vessel Length
Forty-five rabbits were subjected to the patch procedure. Three rabbits at each time point (control, 1-14 days after the patch procedure) were euthanized; the patch-bearing left the common jugular vein and its surrounding tissues were subjected to histological evaluation. The tissues were cut into six 5 mm thick serial slices in an approximately horizontal plane including the cross sections of the arterial patch and the second branch of the left subclavian artery, which were embedded in paraffin. The paraffin-embedded samples were subsequently cut into 4 µm thick sections and stained with Haematoxylin-Eosin and Elastica-Masson. Glass slides were converted into virtual slides using the NDP slide scanner (Nanozoomer 2.0HT) (Hamamatsu Photonics K.K., Shizuoka, Japan), a virtual microscopy system. Measurement of new blood vessel length was carried out using a distant measuring tool of the viewing platform (NDP.View2) (Hamamatsu Photonics K.K., Shizuoka, Japan). Six slices were observed per subject, and 5 to 10 new blood vessels were identified per slice.
We used run_main_code.py to run the whole code and coef_setting.py to adjust parameters.

Conclusions
The process of AVM formation is still unknown. This is because human AVMs are rarely recognized before the onset of symptoms, and their condition can only be confirmed at the time of surgery or autopsy. Furthermore, since there is no animal model that accurately reproduces AVMs, it has been difficult to observe the formation process of AVMs and changes in angiogenic factors over time.
In this study, by observing the AVM animal model newly developed in [10] over time, we confirmed the cyclic variation of VEGF in the early stage of AVM formation. The mathematical model obtained from this study was able to predict the formation of AVMs owing to temporal changes in VEGF. The mathematical simulation reproduced the tangled neovascularization characteristic of AVMs. The simulated elongation process of neovascularization until the formation of arteriovenous fistula suggested that some neovascularization was discarded. This study is expected to shed more light on the early formation process of AVMs, which has not been observed so far.
In addition, in recent years, anti-VEGF monoclonal antibodies have been used clinically in cancer treatment and are widely administered to various carcinomas [28][29][30]. The results of this study suggest that cyclic variation of VEGF plays an important role in the formation of AVMs. Therefore, the administration of anti-VEGF monoclonal antibodies can be considered as a therapeutic agent for AVMs, and it is expected that the effective timing of administration can be predicted mathematically. Acknowledgments: We thank Nobuyuki Takakura from the Research Institute for Microbial Diseases, Osaka University for his assistance in understanding vascular development and maturation through angiogenesis.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A
Before we define our chemotactic velocity v m , we analyze the LEGI (local excitation, global inhibition) model of chemotaxis [31][32][33], and show that three factors play a role: signal gradient, velocity, and temporal change in gradient, all of which are essential as the driving force of cell migration. In the LEGI model, the external (chemoattractant) signal is denoted by S, and the activator and inhibitor derived from signal S are denoted by A and I, respectively. S affects A and I as follows: where k a , k i , γ a , and γ i . The output response of the signal within the cell, denoted as R, is expressed as where F and G are decreasing and increasing functions of R, such as The concrete forms of F and G do not concern the argument below.
In the steady state of (A2), we find that can be described as an increasing function of R. Therefore, we can use the activatorinhibitor ratio Q instead of R to understand the output response of the cell.
Next, we approximate the solution of A(t) and I(t). We define we get local uniformly in t, f ∈ C 1 (R), and also local uniformly in t > 0. This means that A(t) can be approximated by f (t − ) under the assumption γ a 1. That is, The approximation of I(t) can be obtained in a similar manner. The cell has the front and the back: x = ±l. Thus, we have A ± = A(±l, t) and I ± = I(±l, t). With the approximation in (A10), we have dA + dt ∼ γ −1 a k a S t (l, t − γ −1 a ) ∼ γ −1 a k a · γ a S(l, t) − S(l, t − γ −1 a ) ∼ k a S(l, t) − γ a A + (t) (A11) as the approximation of (A1). A − and I pm can be obtained in the same manner. LEGI means that the activator and inhibitor are assumed to be located locally and globally, respectively, on the surface of the cell. Then, A ± and I ± are approximated as A ± ∼ γ −1 a k a S(±l, t − γ −1 a ), where S av (t) = 1 2 (S(l, t) + S(−l, t)).
By simple calculations, the activator-inhibitor ratio, Q = A/I, which dictates the response R, under conditions (A4) and (A12), and with the assumption γ a , γ i 1, given by where Q 0 = (k a γ i )/(k i γ a ). We provide the proof below. Using Equations (A4) and (A12), we have We let (k a γ i )/(k i γ a ) ≡ Q 0 . We may assume the following approximation, using the Taylor expansion. From (A13), we approximate S av (t) ∼ 1 2 (S(l, t) + S(l, t) − 2lS x (l, t)) = S(l, t) − lS x (l, t) (A18) S av ∼ 1 2 (S(−l, t) + 2lS x (−l, t) + S(−l, t)) = S(−l, t) + lS x (−l, t) (A19) for the calculation of x = −l. Applying (A16)-(A18) into (A15) yields Q + (t) ∼ 1 − γ −1 a S t (l,t) S(l,t) 1 − lS x (l,t)−γ −1 i (S t (l,t)+lS xt (l,t)) S(l,t) ∼ Q 0 1 − γ −1 a S t (l, t) S l,t 1 + lS x (l, t) − γ −1 i (S t (l, t) − lS xt (l, t)) S(l, t) ∼ Q 0 1 + −γ −1 a S t (l, t) + lS x (l, t) − γ −1 i (S t (l, t) − lS xt (l, t)) S(l, t) Similarly, For the cell to move in a positive direction, we must have This shows that the chemotaxis is determined by the three factors S x , S t , and S xt , which correspond to the chemoattractant gradient, velocity, and temporal change in the gradient, respectively. The first two factors S x and S t were mentioned in [32]. The additional cell polarity in the direction of S xt inspired us to use this temporal change in the chemoattractant gradient as the chemotactic velocity.