Numerical Approach to a Nonlocal Advection-Reaction-Diffusion Model of Cartilage Pattern Formation

: We propose a numerical approach that combines a radial basis function (RBF) meshless approximation with a ﬁnite difference discretization to solve a nonlinear system of integro-differential equations. The equations are of advection-reaction-diffusion type modeling the formation of pre-cartilage condensations in embryonic chicken limbs. The computational domain is four dimensional in the sense that the cell density depends continuously on two spatial variables as well as two structure variables, namely membrane-bound counterreceptor densities. The biologically proper Dirichlet boundary conditions imposed in the semi-inﬁnite structure variable region is in favor of a meshless method with Gaussian basis functions. Coupled with WENO5 ﬁnite difference spatial discretization and the method of integrating factors, the time integration via method of lines achieves optimal complexity. In addition, the proposed scheme can be extended to similar models with more general boundary conditions. Numerical results are provided to showcase the validity of the scheme.


Introduction
A central question in several areas of biology is how populations of cells interact to form various complex structures. Examples include organogenesis in developmental biology [1,2] or tumor growth in cancer research [3,4]. These phenomena typically span several spatial and temporal scales and involve interactions of many different processes, such as undirected or directed cell motion, cell-cell adhesion, gene expression and the transport of various molecules by diffusion or advection.
Mathematical models have played an important role in the investigation of the underlying cellular and molecular mechanisms of these phenomena. In the framework of PDEs, coupled systems of advection-reaction-diffusion equations are often used. Several modeling approaches for cell-cell adhesion have been proposed in this continuous framework [5][6][7][8], the most widely adopted being that of Armstrong et al. [5] based on a nonlocal adhesion flux term. Most prominently, this has been applied in the area of tumor modeling and cancer invasion, see, e.g., [9][10][11][12][13][14][15]. A similar approach was also used to model the formation of liver cell aggregations in vitro [16].
In this paper, we study an advection-reaction-diffusion system modeling the formation of cell aggregations in cultures of embryonic chicken limbs proposed in [17]. These aggregations, also called pre-cartilage condensations, occur in vivo in the limb bud in early stages of development and lay down the precursors for the skeleton. The model contains nonlocal advection and reaction terms and is effectively four dimensional, with two spatial dimensions and two additional dimensions given by two structure variables on which the cell density depends.
High dimensionality is somewhat typical of models describing the interaction of cell populations. Indeed, it is often necessary to consider not only two or three spatial dimensions, but also additional structure variables for the cell population, such as cell age (modeled together with cell-cell adhesion in [18] and with haptotaxis, e.g., in [19]), or cell size [20]. Taking into account the nonlocal advection terms and the nonlinearity of the reaction terms, as well as the stiffness from highly distinct diffusivities, the numerical simulation of these models is challenging. In fact, due to the computational expense of the model in [17], the original paper only considered numerical solutions in one spatial dimension using a Lax-Friedrichs finite difference scheme. For a thorough computational analysis in the proper biological setting with two spatial dimensions, more sophisticated numerical methods are needed. This is the motivation of this work.
In the rest of this section, we will briefly describe the model from [17], followed by an introduction of the corresponding numerical approach.

Mathematical Model
During early limb development in vertebrates, a template for the skeletal structure is laid out via the aggregation of mesenchymal cells as the first step of chondrogenesis, the formation of cartilage. These cells form tight aggregations called pre-cartilage condensations. At later stages, these condensations differentiate into cartilage, and much later typically into bone. This process can be observed in cell cultures as well. For this, chicken or mouse mesenchymal cells are extracted from early limb buds, mixed and plated as in vitro cultures, so-called high density micromass cultures. In these cultures, cells spontaneously form tight aggregations just like in vivo [21][22][23]. Indeed it is widely accepted that this process recapitulates the process of pre-cartilage condensation formation in vivo [24]. The aggregations in micromass experiments show a variety of morphologies, ranging from spot-like to labyrinth-like patterns, depending on whether the cells were taken from hind or front limbs, as well as various other parameters of the experiments [23]. Very important among these is the initial density of cells with lower cell densities producing spot-like patterns and higher densities producing stripe-like or labyrinth-like patterns [25]. See Figure 1 for some examples from the experimental literature. To elucidate this complex interaction, the following mathematical model was proposed in [17]: with Figure 1. Chondrogenesis patterns in vitro. Image of typical in vitro patterns of pre-cartilage condensations in micromass cultures of mesenchymal cells from chicken or mouse embryos. The dark regions are sites of aggregated cells. Images illustrate different morphologies from spot-like to stripe-like to labyrinth-like patterns. Sources are given with data of initial cell density in the growth medium. (Sources are [25][26][27].) Several mathematical models have been proposed to uncover the mechanism by which pre-cartilage condensations form (see reviews [28][29][30][31][32] and the papers cited therein). Of special interest have been reaction-diffusion type models which give rise to patterns via Turing-type mechanisms, e.g., the recently suggested BSW model [33].
In [17], the authors propose a model of limb chondrogenesis based on experimental work uncovering the role of a network of galectins in in vitro cultures of chicken mesenchyme [34]. Galectins are a class of sugar binding proteins. Galectins can promote or inhibit cell-cell adhesion, a crucial process for the formation of cell aggregations. Galectins are not membrane-bound, but diffuse freely throughout the extracellular matrix. They bind to membrane-bound proteins called counterreceptors. There is a complex interaction between galectins bound to counterreceptors and the production of galectins by cells.
In short, there are two galectins, CG (Chicken Galectin)-1A and CG-8, and two counterreceptors, CR-1 and CR-8. CG-1A may bind to CR-1; the complex induces the production of CG-8 and also crucially mediates cell-cell adhesion. The complex of CG-8 with its counterreceptor CR-8 induces the production of CG-1A, so that the two galectins effectively give rise to a positive feedback loop. However, CG-8 can also bind to the counterreceptor CR-1, thus 'blocking' these for CG-1A. So somewhat paradoxically, an increase of CG-8 leads to fewer condensations while an increase of CG-1A, enhancing cell-cell adhesion, leads to more and tighter condensations. (See [34] for further details about the biochemical interactions.) To elucidate this complex interaction, the following mathematical model was proposed in [17]: with Here R(t, x, T 1 , T 8 ) denotes the cell density as a function of time, position and membrane bound counterreceptor densities T 1 and T 8 (for CR-1 and CR-8, respectively). The concentrations of the galectins CG-1A and CG-8 are given by c u 1 (t, x) and c u 8 (t, x). (See Table 1 for a short list of the variables used.) We assume that x ∈ Ω, which is a finite spatial region with piecewise smooth boundary in the plane. The parameter region is denoted by Ω T = [0, ∞) × [0, ∞) under the assumption that T 1 ≥ 0 and T 8 ≥ 0. Differentiation with respect to the spatial variable x is denoted by ∇ x .

Variable
Meaning density of counterreceptor CR-1 (membrane-bound) T 8 density of counterreceptor CR-8 (membrane-bound) R(t, x, T 1 , T 8 ) cell density as a function of time, space, counterreceptor densities c u The model is of reaction-diffusion-advection type, where the spatial advection term models cell-cell adhesion; we refer to this as a reaction-diffusion-advection system. The reaction terms in (2) and (3) are nonlocal in order to model the dependence of galectin production on the counterreceptor density on the cell membranes. The rates (6) and (7) model the production of counterreceptors. Crucially, cell-cell adhesion is modeled via a nonlocal flux term in (1), following the cell-cell adhesion model by Armstrong et al. [5]. The essential idea is that the effective cell flux term K(R) induced by cell-cell adhesion is the result of a weighted averaging process. Cells at position x move in the direction of greatest cell density within a disk D r 0 (x) of radius r 0 , the sensing radius, mediated by the density of the CG-1A/CR-1 complex, denoted by c 1 . Here we use a logistic form for σ in the expression for the adhesion flux, i.e., The parameter R m is a maximum density. In (8), D r 0 (x) denotes a disk of radius r 0 centered at x in two spatial dimensions. We use periodic boundary conditions on the spatial domain ∂Ω and homogeneous Dirichlet on ∂Ω T throughout this work. The use of periodic boundary conditions is motivated by the fact that we are interested in modeling the behavior of the cell population away from the boundary, a common approach for models of cell populations to exclude boundary effects, both in PDE models (e.g., [5,35,36]) and cellular automata models such as the Cellular Potts Model [37], see, e.g., [38]. The proposed numerical methods can be extended to other biologically meaningful boundary conditions, most prominently no-flux boundary conditions on the spatial domain.
The model was analyzed numerically in one spatial dimension in [17]. It was shown to be able to produce repeating patterns in the cell density R and to explain several experimental results. In this work, we propose a numerical method for two spatial dimensions and present related numerical results on the behavior of the system, notably that the model can produce a range of patterns which qualitatively agree well with the experimental patterns, see   Note that below a critical curve in R m α K -plane, the system does not produce patterns, and that above this curve, labyrinth-type, stripe-like, or spot-like patterns may arise, similar to the experimental data in Figure 1. Other parameters used were as in Table 5.
A variant of the model was also proposed in [17], called the "reduced system". This was derived from the model (1)-(8) based on the assumption of fast counterreceptor production, i.e., that the timescale of production of galectins is much longer than that of counterreceptors. This effectively eliminates the dependence of the cell density R on T 1 and T 8 , so that R is only a function of time and space; below we denote by R = R(t, x) the total cell density. The reduced system is now Here K(R) can be computed using (8) which yields where a logistic dependence of the flux on the cell density as in (9) is used, namely The reduced system (10)- (12) is defined on the three dimensional space Ω t × Ω, where Ω t ⊆ R is a time interval. This is in contrast to the original Equations (1)- (8) , which is defined on the five dimensional space Ω t × Ω × Ω T . This difference in dimensionality makes the numerical treatment of the reduced system much easier. In addition, it has a spatially homogeneous steady state solution, and thus standard linearization techniques are applicable to investigate its stability and the characteristic wavelength of patterns it gives rise to.
Biologically, it is unclear whether the assumption of fast counterreceptor production, underlying the derivation of the reduced system, is realistic since data on the relative time scales is not available [17]. Numerical results in one spatial dimensions indicate that the two systems (10)- (12) and (1)-(8) display qualitative similar behavior and that the characteristic wavelength of the reduced system (10)- (12) can be used to predict that of the full system (1)- (8). However, there are indications that the reduced system cannot adequately capture the subtle inhibitory role the galection CG-8 has on condensation patterns: When it is added in vitro, the wavelength of the resulting pattern increases. This is the case for certain parameters of the full model, but not of the reduced model [17].

Numerical Approach
The proposed numerical approach is based on a variable decoupling treatment to (1). For fixed x, R(t, x, T 1 , T 8 ) is governed by an advection equation on Ω T with Dirichlet boundary conditions. Whereas for fixed (T 1 , T 8 ), R(t, x, T 1 , T 8 ) is governed by an advection-diffusion equation on Ω with periodic boundary conditions. These two processes can be handled differently with a mesh-free scheme and a mesh-based scheme respectively. We briefly describe the discretization strategies in Ω T and Ω.
At each time step, the major numerical complexity comes from evaluating the four-variable function R(t, x, T 1 , T 8 ).
Since Ω T is unbounded, a radial basis function (RBF) meshless approximation is combined with the regular spatial finite difference scheme to complete the 4D discretization of R. That is, R is expanded as a linear combination of RBFs for fixed x while being discretized by regular finite difference scheme for fixed T 1 and T 8 . Hence all the double integrals over Ω T as well as the partial derivatives in T 1 and T 8 can be evaluated exactly. K(R) can then be evaluated using Simpson's rule.
The divergence ∇ x · (RK(R)) is computed using a 2D finite difference fifth order weighted essentially non-oscillatory scheme, the popular WENO5 [39], so that a non-oscillatory and sharp spatial profile of R maintains during the pattern formation.
The time integration is based on the method of lines. In the time iteration, a large stiffness ratio exists due to the highly distinctive diffusivities in the cell density R and the concentrations of the galectins CG-1A and CG-8, c u 1 and c u 8 . The method of integrating factors (IF) [9,40,41] is applied to transform the original system of PDEs into ODEs which can be solved by a proper ODE solver. With periodic boundary conditions, the evaluation of matrix exponents involved in the integrating factor method can be accomplished via fast Fourier transform (FFT). This makes IF an optimal choice in the sense that diffusion operators can be integrated exactly and effortlessly.

Outline of Paper
This paper is organized as follows. In Section 2, we provide an overview of the related numerical schemes, including the RBF meshless method, WENO5 scheme, the integrating factor method and a third order time integration method. The numerical algorithm and the corresponding complexity analysis will also be provided. Section 3 is devoted to the numerical results that validate the proposed numerical approach. Finally, concluding remarks and future investigations will be given in Section 4.

Numerical Implementation
To make this paper self-contained, we start with a brief overview of the numerical schemes involved.

Radial Basis Function (RBF) Expansion
For fixed t and x, R as a function of T 1 and T 8 is defined over the semi-infinite region Ω T = [0, ∞) × [0, ∞). A RBF meshless approximation will be much more cost-effective than the regular finite difference discretization. To make the model biologically meaningful, homogeneous Dirichlet conditions are imposed over ∂Ω T which suggest the choice of radio basis functions in the form of They satisfy the far end boundary conditions automatically. Here we use T = (T 1 , T 8 ) and T j = (T j 1 , T j 8 ) to represent points in the T 1 T 8 -plane. In particular, T j , ∀j = 1, · · · , M 0 , are the M 0 source points used to define the M 0 radio basis functions. α > 0 is a scaling factor termed shape parameter that can be adjusted to enhance the numerical accuracy, which will become clear later. · stands for the Euclidean norm.
At each iteration step and fixed spatial location, the T 1 and T 8 dependence of R is denoted by where q j , ∀j = 1, · · · , M 0 , are the coefficients of the RBF expansion to be determined.
Once the grid values of R have been calculated, we determine the coefficients via interpolations. Let M I be the number of interior sample points and M − M I be the number of boundary sample points. They are also termed collocation points. Fitting the values of R at the interior pointsT i ∈ (0, ∞) × (0, ∞), ∀i = 1, · · · , M I and the boundary pointsT i on the edges T 1 = 0 and T 8 = 0, ∀i = M I + 1, · · · , M, yields the following M × M 0 linear system for the coefficients q j s, whereR is an M I dimensional vector storing the numerical values of R at the M I interior sampling locations, and 0 is the (M − M I ) dimensional zero vector reflecting the homogeneous Dirichlet conditions. The coefficient matrix of the above linear system is dense but of a small scale. For α > 0 and M = M 0 , it is positive-definite [42] with the Gaussian basis functions, which guarantees the uniqueness of the solution. Desired condition number of the coefficient matrix can also be achieved by properly adjusting the scaling factor α, the source and collocation points, T j andT j [42]. Least squares approximation can be used to solve the M × M 0 linear system when M > M 0 , which is preferred in this work.

Finite Difference WENO5 Discretization
The spatial divergence ∇ x · (RK(R)) is determined using WENO5 [39]. For simplification of notation in the following illustration, we temporarily drop the c u 1 dependence and denote This is valid since R will be updated for fixed c u 1 and c u 8 in the numerical simulation. Consider the spatial discretization for any fixed (T 1 , T 8 ). Once the numerical values of R are computed over a uniform 2D mesh in the spatial domain, the grid values of the nonlinear flux term F(R) can be generated via the RBF approximation and a high-order numerical integration scheme such as Simpson's rule. Then the 2D finite difference WENO5 reconstruction based on the grid values F(R ij ) can be implemented to output the numerical fluxesF 1 Define the smoothness indicators as The convex combination is applied to determine the numerical flux at the cell center. Here The constants involved are set as to determine the spatial discretization of the nonlinear flux term ∆y .
Please note that obtaining α 1 and α 2 analytically is quite complicated due to the integral form of the nonlinear flux term. Numerical approximation of them is sufficient to maintain the monotonicity.

Time Integration
We first introduce the following simplified notation for the full model.

Method of Integrating Factors
Using the method of integrating factors, the system of PDEs (15)- (17) are transformed into the following ODEs where e −tµ∇ 2 denotes the operator exponential involving the diffusion operator with diffusivity µ.

Time Integration Method
Equations (18)- (20) are all in the general form of which can be solved using various time integrators [43][44][45][46]. We use a third-order explicit Runge-Kutta method (SEIF3) [47]. Denote the numerical approximation of u at time step n by U n and the spatial discretization of µ∇ 2 by C. Then SEIF3 reads With periodic boundary conditions on ∂Ω, the evaluation of e β∆tC U for any matrix β∆tC can be achieved using FFT. In case of non-periodic boundary conditions, variations of the presented IF coupled with FFT can be applied, such as compact integration methods via Krylov subspace approximations [43][44][45][46]. The main idea is to approximate the resulting vector e β∆tC U by a vector in the span of some properly chosen basis that depends on the matrix C and imposed boundary conditions. The FFT approach based on Fourier series presented here can be viewed as a special case of the general Krylov expansion.

Numerical Algorithm and Complexity
The detailed numerical algorithm is provided as follows.
• Let M 0 be the number of source points in the T 1 T 8 -plane. Initialize M 0 different R's over a 2D spatial mesh of N × N grids, associated with the source points. • Initialize c u 1 and c u 8 over the same 2D spatial mesh of N × N grids. • Compute exactly the partial derivatives and integrals of radio basis functions to be used later, including • Perform the following at each time step up to the desired t = T.

1.
Determine the coefficients q j 's in the RBF expansion for all the R's.

4.
Determine G 1 (R, c u 1 , c u 8 ) and G 8 (R, c u 1 , c u 8 ) via linear combinations of the known information from the radio basis functions.

5.
Update c u 1 and c u 8 explicitly using SEIF3.
Steps 1 and 3 are most expensive.
Step 1 requires solving N 2 linear systems of size M × M 0 (M ≥ M 0 ) via least squares with QR factorization. The complexity is O(M 3 0 + N 2 M 2 0 ). For each (T 1 , T 8 ) block in Step 3, the worst-case complexity of FFT is O (N 2 ln N), hence the total complexity is O(M 0 N 2 ln N) for N M 0 . The efficiency of the RBF meshless approach is due to its rapid convergence which requires rather small value of M 0 .

Numerical Results
In this section, we present some numerical results for both the reduced and full systems. The biological parameters are as shown in Table 5. (See also references [17,48].) In all the following simulation, the spatial domain is set to be Ω = [0, 1] × [0, 1] so the rectangular spatial meshing is a routine process. Nevertheless the domain of computation is not confined to rectangles. For example, a circular region can be transformed into a rectangular one in polar coordinates hence all the numerical discretizations will still be valid with a change of variables. In case of general irregular domains, finite element WENO or the recently developed RBF-WENO can be applied instead of the finite difference WENO used in our work. This is one of our follow-up works in process.

Stability of Time Integration Method
The reduced model is solved to confirm the stability of the proposed time integration method coupled with WENO5 spatial discretization.
Since c u 1 and c u 8 diffuse much more rapidly than R within the time scale of simulation, explicit time iteration for the reaction terms is sufficient to capture the mild changes of them without stringent time stepping restrictions. Initialize R, c u 1 and c u 8 close to their spatially constant steady states. Set the uniform spatial grids with mesh size ∆x = 1 2 7 . Time step ∆t = 0.3∆x guarantees the stability of SEIF3 for α K = 60, 80 and 100, as shown in Figure 3 where the time evolutions of the total variation (TV) of R for these parameter choices are plotted. Denote the ij-th spatial grid by (x i , y j ). The total variation is computed numerically via Starting from the beginning, TV of R increases rapidly due to the dominance of cell-cell adhesion over diffusion, which leads to the gradual formation of patterns until TV reaches a peak. Then the adhesive force is well balanced by the diffusion and TV maintains afterwards except for the slight fluctuations due to the shape adjustment of patterns at the final stage. Clearly pattern formation occurs more rapidly with higher α K to which more detailed biological explanation will be provided in Section 3.4. The conservation of spatial mean of R throughout the time iteration can also be confirmed by our computation. As an example, the evolution of spatial mean of R with α K = 60 is shown in Figure 4. . Plots of total variation of the cell density R versus time (0 ≤ t ≤ 2) in the reduced system (10)- (12). Cases with α K = 60, 80 and 100 are compared. Other related parameters are given in Table 5.  Table 5.

Convergence of RBF Approximation
Recall that R(t, x, T1, T8) is governed by an advection equation with homogeneous Dirichlet boundary conditions at each fixed spatial location x; see (1). To study the dynamics of R in T 1 T 8 -space, let us temporarily drop x. Advection in T 1 T 8 -space is determined by the velocity field ( γ(c u 1 , c u 8 , T 1 ), δ(c u 8 , T 8 )) given in (6) and (7), respectively. Keeping x and c u 1 (t, x), c u 8 (t, x) fixed, let us write γ(c u 1 , c u 8 , T 1 ) = γ(T 1 ) and δ(c u 8 , T 8 ) = δ(T 8 ). A plot of the velocity field (γ, δ) in the T 1 T 8 -plane is given in Figure 5. There is a unique zero (T 1 , T 8 ) with T 1 , T 8 > 0, and every flow line with positive initial T 1 converges to (T 1 , T 8 ). This means that R is advected towards this zero and, again fixing the position x and the concentrations c u 1 (t, x), c u 8 (t, x), R becomes concentrated at (T 1 , T 8 ). In fact, under simplifying assumptions, it is possible to show that R(t, T 1 , T 8 ) converges to a Dirac delta distribution centered at (T 1 , T 8 ) [17]. This is the assumption in the case of fast counterreceptor dynamics (i.e., fast dynamics in T 1 T 8 -space), which yields the reduced model (10)- (12). Thus, the full model is relevant for the case that the flux term in T 1 and T 8 takes on values for which the variance of the distribution R(t, T 1 , T 8 ) is large enough that it cannot be reasonably approximated by a delta-distribution, i.e., that the decay from the peak is not too rapid to lose the smoothness of R(t, T 1 , T 8 ). On the other hand, the decay is rapid enough so that the support of R(t, T 1 , T 8 ) can be approximated by a compact computational region, which is set as At each fixed spatial location x, we initialize R with an Gaussian distribution This assumption of a joint Gaussian distribution for the membrane-bound counterreceptor densities T 1 and T 8 is biologically plausible. Indeed, since the density of counterreceptors on an individual cell may be modeled as the cumulative outcome of many, partly uncorrelated small-scale events, an overall normal distribution is a reasonable assumption by an appeal to the central limit theorem; this standard argument is well laid out in e.g., [49]. Empirical data supports this normality assumption, see, e.g., Figure 5 in [50] showing the distribution of PTK7 receptor density for human HeLa cells. The above heuristic arguments mean that a series of deformed Gaussian distributions centered at various (T 1 , T 8 ) locations, R(t, T 1 , T 8 ), will be generated through the advection process and then determined via a RBF approximation. That is R(t, T 1 , T 8 ) will be written as a linear combination of the Gaussian basis φ j introduced in Section 2.1 for t > 0. The accuracy of the RBF approximation to R(t, T 1 , T 8 ) depends on the variance and center location of R(t, T 1 , T 8 ), the source and sample points as well as the shape parameter α in the Gaussian basis. Recall that the source points are the data points used to define the RBFs and the sample points are the data points used for interpolation. They are set to be the lattice grids as shown in Figure 6 to fit the rectangular domain of computation. To properly choose the shape parameter α, we shall study the effectiveness of the RBF approximation to various Gaussian functions since they behave similarly as R(t, T 1 , T 8 ) based on the above discussions. Here are the constructive observations from numerical simulation.

•
Decreasing α flattens the Gaussian basis functions but significantly increases the condition number of the coefficient matrix in the linear system (14), some examples are given in Table 2.  • A Gaussian function with shape parameter α 0 is more accurately approximated by the Gaussian basis functions with shape parameter α closest to α 0 , as shown in Table 3. • The Gaussian basis approximation to a flatter Gaussian function achieves higher order of accuracy with matching Gaussian basis functions (i.e., with α closest to α 0 ), as shown in Table 3. • A Gaussian function centered close to the boundary is less accurately approximated by any Gaussian basis, as shown in Table 4. This effect is most obvious for the flattest Gaussian function (left most column in Table 4) since it does not decay fast enough to smoothly match the homogeneous Dirichlet condition at boundaries T 1 = 0 and T 8 = 0. When the variance of a Gaussian function decreases, it decays fast enough so that the boundary disagreement is negligible compared to the interpolating error caused by the non-flatness of the Gaussian function. This explains why the minimum error does not have to occur exactly at the center of [0.8] × [0, 8], referring to the right most column in Table 4.
Numerical results of R(t, T 1 , T 8 ) will be provided in Section 3.4 in which we can observe that the mean peak location of R(t, T 1 , T 8 ) stays reasonably away from the boundaries. As a result, setting α = √ 10 in all the numerical simulation yields desired conditioning and approximation accuracy.   Table 3. RBF approximation errors in interpolating different Gaussian functions with shape parameters α 0 = 1, √ 10 and √ 60 using Gaussian basis functions with shape parameters α = 1, 2, · · · , 10, all centered at (4.5, 4.5). The error is defined as the absolute difference between the mean values of the original Gaussian function and its RBF approximation. Regardless of α 0 , the error is minimized when α is closest to α 0 . In the left most column, the minimum error is of smallest scale when the flattest Gaussian function among the three is approximated. However it is most sensitive to the change of α due to the ill-conditioning of the coefficient matrix.  Table 4. RBF approximation errors in interpolating different Gaussian functions with shape parameters α 0 = 1, 2, 3 and 4 and centered at integer locations (k, k), for 1 ≤ k ≤ 8, using Gaussian basis functions with shape parameters α = √ 10. The error is defined as the absolute difference between the mean values of the original Gaussian function and its RBF approximation. Regardless of α 0 , larger approximation error is observed when interpolating a Gaussian function centered close to the boundary.

Numerical Examples for the Full System
Starting from initial conditions close to a spatially constant distribution of the cell density, the system forms characteristic stripe-like or labyrinthine patterns similar to patterns observed in cultures of embryonic chick mesenchymal cells in vitro; see Figure 1 for photos of corresponding experimentally obtained patterns. Figure 7 illustrates the process of pattern formation in the full model (1)-(3) as labyrinthine patterns in the total cell density emerge. These patterns are quasi steady state patterns in the sense that after their rapid establishment, they change little over long time scales. Please note that the overall concentrations of CG-1A and CG-8, c u 1 and c u 8 , are increasing throughout the simulation. This is in accordance with experiments [17,34].  Table 5, in particular f = 0.8.

Mechanism of Pattern Formation
The term responsible for the formation of patterns in the total cell density ∞ 0 ∞ 0 R(t, x, y, T 1 , T 8 )dT 1 dT 8 is the cell-cell adhesion term K(R) in (8). This term describes an effective motion up the gradients of the cell density R, but modulated by the CR-1 counterreceptor density T 1 , so that the larger T 1 , the larger the flux. This flux is balanced by the diffusive flux in the quasi steady state. An illustration of the adhesive flux is given in Figure 8.  Figure 7 is shown for the final time t = 1.5. A plot of the adhesion flux K(t, x, T 1 , T 8 ) for T 1 = 2.0 and T 8 = 3.5 in red is displayed along with a density plot of the total cell density. Right: Detail of the contour plot to K(t, x, T 1 , T 8 ). Please note that adhesion generates an effective velocity up the gradients of total cell density, which is responsible for the aggregation of cells. In the quasi steady state, this flux is balanced by the diffusive flux (not shown).
Because the cell adhesion flux term K(R) increases with counterreceptor density T 1 , the cell aggregations (areas of high cell density) are primarily made up of cells with large density T 1 . This can be seen in the plots of the mean of T 1 in Figure 7. Regions of larger mean T 1 also correspond to regions of higher total cell density. This observation is confirmed when plotting the cell density against the mean counterrecptor density T 1 in Figure 9. Initially (time t = 0), there is no relation between these two variables since the cells are assumed to be completely mixed in the initial data. At the end of the simulation (t = 1.0), there is a clear correlation between cell density and mean T 1 . Figure 9. Plots of total cell density R tot = R(t, x, T 1 , T 8 )dT 1 dT 8 versus the mean of the counterreceptor T 1 , given by T 1 R(t, x, T 1 , T 8 )dT 1 dT 8 /R tot for initial time t = 0 (left) and final time t = 1 (right). Please note that there is initially no correlation between the cell density and the counterreceptor (by construction of the initial conditions), but in the final pattern, regions of high cell density are made up of cells with higher CR-1 counterreceptor density T 1 . (The data is the same as in Figure 7.) The counterreceptor CR − 1 not only mediates cell-cell adhesion, but it is also indirectly involved in the production of CG-1A, so that this numerical result is consistent with the experimental finding that production of CG-1A is elevated inside condensations [17,34].

Comparison of Solutions to the Full and Reduced Models
Since the cell density R(t, x) in the reduced system (10)- (12) does not depend on the counterreceptor densities T 1 and T 8 , the related simulation is only two-dimensional compared to the effectively four-dimensional full system. This is certainly a computational advantage over the full system (1)-(3). However, its derivation rests on the additional modeling hypothesis that production of counterreceptors happens on a faster time scale than production of galectins, an assumption whose validity is unknown and indeed difficult to test experimentally [17]. While the system still produces patterns in the cell densities, and thus can be seen as a simple, but useful model of chondrogenesis, there are certain aspects of the full model that are lost. For instance, the parameter f describes the affinity of CG-8 to bind to the counterreceptor CR-1 relative to the counterreceptor CR-8 and thus is vital to an understanding of the role of CG-8. By binding to CR-1, CG-8 is effectively in competition with CG-1A; thus f describes the 'competitiveness' of CG-8. In the reduced model, f does not appear because in the steady-state assumption of fast equilibration of CR-1, the density c 1 of the complex of CG-1A and CR-1 is independent of f . Thus, the reduced model is inadequate to investigate the subtle role of CG-8. As an illustration, we show some numerical results in Figure 10. All three simulations are corresponding to the same parameter set, with the exception of the parameter f . In the case of f = 0.8, or high competitiveness of CG-8, a clear pattern in the total cell density emerges. Areas of high cell density also correspond to areas of high average density of CR-1 (see Figure 7, column 2). In other words, the condensations are formed by cells with higher average CR-1 densities on their membranes.
In the case of f = 0.2, or 'low' competitiveness of CG-8, the cell pattern is even more pronounced with very clear boundaries between high and low cell densities . This is because cell-cell adhesion is mediated by the complex of CG-1A and CR-1, which is now less affected by the competition of CG-8, leading to stronger cell-cell adhesion. Finally, in the reduced model, the dependency on the parameter f is lost through the steady-state assumption for CR-1. Now all cells have essentially the same density of CR-1. The pattern is much fainter, pointing to the crucial effect of different densities of CR-1, i.e., the dependence of the cell density on the variable T 1 . These examples illustrate that the full model enables a much more in-depth analysis of the interplay of CG-1A and CG-8 and its significance for the formation of condensations.   (12). (The reduced model is independent of the parameter f , see text above.) The initial conditions for the two full-model simulations were identical. Note the different scales for the cell densities. Higher f means higher 'competitiveness' of CG-8, leading to reduced cell-cell adhesion and fainter condensations. The reduced model also shows condensations, but much fainter, pointing to the importance of variability in the dependence of the cell density R(t, x, T 1 , T 8 ) on the variable T 1 in the full model.

Conclusions and Future Investigations
We present a RBF meshless approximation coupled with the mesh-based finite difference discretization to relax the computational complexity caused by the high dimensionality in a nonlinear advection-reaction-diffusion system. The span of rather small sets of Gaussian basis functions can provide desired accuracy in the numerical approximation of the structure variable dependence of the cell condensation to be captured. This high efficiency of the meshfree approach significantly reduces the computational cost of a 4-dimensional system to one that is comparable to a 2-dimensional system. We presented some sample simulations that illustrate how this system gives rise to spatial patterns in the cell density resembling in vitro condensations. Additional insights into the mechanism of pattern formation can be gained by investigating not only the final pattern, but also the temporal development of the process, and also by simulating specific experiments [34,48]. We will present these results in a separate forthcoming publication.
The RBF meshfree approximation not only couples well with the finite difference spatial discretization presented in this work, but also will fit in other types of spatial discretization, such as finite element or finite volume methods, depending on the properties of spatial domain as well as the prescribed boundary conditions. Furthermore, a complete meshfree scheme in both structural and spatial domains will be one of our future attempts, allowing applications to extensions of the chondrogenesis model at hand as well as other biological models. This is especially relevant for modeling in vivo chondrogenesis, which takes place within the growing limb bud, thus necessitating growing domains with curved boundaries [51,52].
Another numerical extension of this work is to optimize the RBF approximation error through the development of adaptivity strategies for RBFs via dynamically adjusting the source points, the sample points and the shape parameters. This will allow cost-effective meshfree approximations of non-smooth or rapidly changing biological quantities involved in various complex dynamical systems.
Author Contributions: Conceptualization, T.G. and J.Z.; Investigation, T.G. and J.Z.; Software, T.G. and J.Z.; Writing-original draft, T.G. and J.Z.; Writing-review & editing, T.G. and J.Z. All authors have read and agreed to the published version of the manuscript.