A mathematical modeling technique to understand the role of decoy receptors in ligand-receptor interaction

The ligand-receptor interaction is fundamental to many cellular processes in eukaryotic cells such as cell migration, proliferation, adhesion, signaling and so on. Cell migration is a central process in the development of organisms. Receptor induced chemo-tactic sensitivity plays an important role in cell migration. However, recently some receptors identified as decoy receptors, obstruct some mechanisms of certain regular receptors. DcR3 is one such important decoy receptor, generally found in glioma cell, RCC cell and many various malignant cells which obstruct some mechanism including apoptosis cell-signaling, cell inflammation, cell migration. The goal of our work is to mathematically formulate ligand-receptor interaction induced cell migration in the presence of decoy receptors. We develop here a novel mathematical model, consisting of four coupled partial differential equations which predict the movement of glioma cells due to the reaction-kinetic mechanism between regular receptors CD95, its ligand CD95L and decoy receptors DcR3 as obtained in experimental results. The aim is to measure the number of cells in the chamber’s filter for different concentrations of ligand in presence of decoy receptors and compute the distance travelled by the cells inside filter due to cell migration. Using experimental results, we validate our model which suggests that the concentration of ligands plays an important role in cell migration.

where R b is the bound receptor. The bound receptor R b modulates pseudopod extension 4,5 which is critical in sensing the direction of the movement.
Decoy receptors are some special type of receptors that acts as an inhibitor in cell signalling induced by normal receptor. IL-1R2 was identified as the first decoy receptor in the '90s 14 . Most of the mechanisms of the decoy receptors are yet to be known in the view of cellular biology 15,16 . The subfamily of decoy receptor includes DcR3 17 , D6 18 , Duffy Antigen Receptor for Chemokines (DARC) 19 , etc. Decoy receptor 3 (DcR3) 17,20 , a well-known member of tumor necrosis factor receptor superfamily (TNFRSF), is capable of neutralizing the biological effect of other three tumor necrosis factor superfamily (TNFSF) members, namely Fas ligand (FasL/ CD95L) 21 , LIGHT 22 and TL1A 23 . Based on its neutralizing effects on CD95L, LIGHT, and TL1A, DcR3 can be defined as an immunomodulator 17,20 , since CD95L and LIGHT induce apoptosis and inflammation, and TL1A is a T cell costimulator. DcR3's overexpression was found in various malignant cells such as RCC cell 24 , Gliomas cell 25 , Glioblasomas cell 26 . Experimental findings by Weissinger et al. 25 and Roth et al. 8 suggest that DcR3 binds with CD95L and neutralizes CD95L chemo-tactic activity. As a result, cell migrated area and cell motility were decreased in the presence of DcR3. The decoy receptor's interaction with ligand was captured in a mathematical model in 27 , but the chemo-repellent activity of the bound decoy receptors does not influence the cell migration in that model nor the application of the model is explored well in 27 . These complicated facets of cell motility are exclusively investigated in this framework of partial differential equations. For instance, we measure the migrated cell density in the Boyden chamber's filter for different concentrations of ligand and compute the distances travelled by the cells inside filter. In particular, we investigate mathematically how the cell migration is influenced by the interaction of normal receptors CD95L, decoy receptors DcR3 and ligand CD95L via a mathematical model.
Chemo-taxis and Brownian motion are used to model chemically modulated cell motility following the pioneering work of Kellar & Segal 28 . In 28 , the authors have considered cell-flux at spatial point x and at an instant of time t ( J n (x, t) ) which is dependent on cell density n(x, t) and chemical concentration c(x, t) through the following equation where D n represents random movement of cell and χ n (c) represents the chemo-tactic sensitivity, which is chemoattractant if χ n (c) > 0 for all c and chemo-repellent if χ n (c) < 0 for all c. Since chemo-attractants or chemo-repellents activate cell's surface receptors, it is interesting to explore how the receptor-ligands interaction affects the motility of cells. Sherratt et al. 13 have proposed a three compartment model of cell n, ligand c and bound receptors ρ with some modification of cellflux of (2) to account the receptor-ligand mechanism. Despite the importance of these research and the crucial outcome they provided, their work did not account for cell migration due to decoy receptors. In this study, we have incorporated decoy receptors as chemo-repellent and chemo-attractant to elucidate the effect of these receptors in cell migration. We assume that extracellular complex produced through interaction of receptor and ligand attracts the cells whereas the ligand-decoy complex repeals. These two reaction occurs simultaneously through multiple receptor state in Boyden chamber. The Boyden chamber 3,29 has become a popular tool to estimate the coefficients of random cell walks and chemo-tactic behaviors for various cell types by fitting the experiments with theoretical computations in vitro system. It is classically constructed with two wells separated by a porous filter, which acts as a physical barrier and active migration is the only way to overcome it. A chemo-attractant solution is placed in the lower well while a cell suspension is placed in the upper well. Cells are allowed to migrate inside the filter in response to the chemical distribution for a given time interval. Cell motility is measured from the distance that cell wave fronts reach inside the filter during the incubation time. Sherratt et al. 13 have shown the movement of cells in filter affected by the normal receptor in Boyden chamber. Chen et al. 30 studied cell migration due to receptor-ligand interaction with cell sedimentation in Boyden chamber. Zigmond and Dunn chambers are two other chemotaxis bridge assays which studies cell motility of polymorphonuclear (PMN) leukocytes and fibroblasts respectively. The Zigmond chamber with efficient optical properties generates shallow gradients so that cells could respond to only 1% changes in concentration 31 . The Dunn chamber studies the migration of fibroblasts whose migration rate is 0.42 − 1.25 µm/min and generates a linear gradient within an hour of its initial setting with a gradient half-life of 10 to 30 hours. The Insall chamber a refinement of Dunn chamber, consists of two-well bridge design assay with defined directions of gradient and two different gradient steepness in the same assay. It measures melanoma chemotaxis 31 .
In this work, we mathematically formulate the effect of cell migration in presence of decoy receptors and normal receptors. To the best of our knowledge, no mathematical model has considered cell migration caused by decoy receptors and studied chemo-repellent activity of the bound decoy receptors. We construct a novel mathematical model to study and validate the movement of glioma cells in presence of normal receptors CD95, ligand CD95L and decoy receptors DcR3 mathematically. The novelty of our mathematical model is the inclusion of decoy receptors in cell migration. The goals of our work are to: (i) propose a new mathematical model to understand the role and dynamics of decoy receptor DcR3 in cell migration based upon the classic Keller-Segel www.nature.com/scientificreports/ model 28 , (ii) study the system dynamics in presence and absence of decoy receptors DcR3, (iii) measure the migrated cell density in the Boyden chamber's filter for different concentrations of ligand in response to chemotactic sensitivity of normal and decoy receptors, (iv) compute the distances travelled by the cells inside filter due to cell migration. Here, we consider both decoy receptors as well as normal receptors present for cell migration. We fill the upper well and lower well with cells and ligands respectively and study the dynamics of the system in filter compartment.

Theoretical study of receptor-ligand modeling
Governing equations. We consider the normal receptor-ligand and decoy receptor-ligand interaction to occur simultaneously as per multiple receptor states 2 . Assume that the normal receptors R f and decoy receptors D f bind with the ligands at a rate k a 1 and k a 2 and form the bound normal receptors R b and bound decoy receptors D b respectively. Let the degradation rate of R b and D b are k d 1 and k d 2 respectively. The bound normal receptor R b ( D b ) turn into internalization components at a rate k i 1 ( k i 2 ). The flowchart of the normal and decoy receptor-ligand interaction is shown in Fig. 1a. Using the law of mass action, the pathophysiology is modeled by the following system of Eq. (3) subject to non-negative initial conditions. Note that the net mass-balance of system (3) is negative due to the internalization of bound receptors and decoy receptors 2,13,32 . The fate of bound receptors can vary depending on several factors, including the receptor type, the nature of the ligand, and the number of times they have interacted. In some cases, the internalized receptors can be recycled back to the cell surface as free receptors with the ligand, allowing for continued signaling. However, in other cases, the bound receptors may be targeted for degradation in lysosomes through internalized endocytosis 2,30 , which can lead to decreased signaling and modulation of cellular responses.
To study chemo-tactic sensitivity of the normal receptors and decoy receptors, we spatially extend system (3) by considering the Boyden chamber as spatial domain. The Boyden chamber (see Fig. 1b) consists of three compartments namely, upper well, lower well and filter. We fill the upper well and lower well with cells and ligands respectively. The movement of the cell depends on the interaction of ligands and normal receptors/ decoy receptors which are situated at the surface of the cell. Neglecting the radial variation, we consider the onedimensional geometry of Boyden chamber as the domain for our spatio-temporal model. Let the thickness of lower well, filter and upper well be L l , L f and L u respectively. We consider the upper surface of the filter as x = 0 . Chemical c is defined on −L u < x < L l + L f as it can diffuse throughout the chamber. We are interested to understand the movement of cells in the filter only in presence of decoy and cell receptors. There is an inflow and outflow of cells in the upper and lower surface of filter respectively. Chemicals can diffuse throughout the chamber with no-flux boundary conditions at both end of the chamber. The bound receptors and bound decoy receptors on the cell surface act as chemo-attractant and chemo-repellent. Therefore the cell-flux is modelled by www.nature.com/scientificreports/ n are the number of normal receptors and bound receptors per cell respectively. Parameters χ and µ are the chemo-tactic attractant and repellent coefficient respectively. The spatio-temporal behavior of the system in 0 < x < L f is given by All the variable used in the system (4) is mentioned Table 1. In other compartments, we assume that chemical c flows through diffusion only. c(x, t) is following the reaction-diffusion equation We also need the continuity of chemical flux ( ∂c ∂x ), at the both end of filter x = 0 and x = L f . Let the diameter of the cell be L, then inflow and outflow cell-flux can be modelled in the boundary condition as In the absence of ligand, a certain number of free receptors are typically present on the cell surface. However, when the ligand is present and binds to its receptors on the cell surface, it can trigger a signaling cascade that ultimately results in the upregulation of additional receptors 30,33 . This upregulation can occur through several mechanisms, including increased transcription and translation of receptor genes, decreased degradation of existing receptors, and increased recycling of internalized receptors back to the cell surface. To incorporate such upregulation mechanism, we assume that total number of free receptors and bound receptors per cell is increasing with rate f (ρ r ) = Ŵ r (1 + βρ r ), where Ŵ r is the initial number of moles of free receptors on a cell and β is the increasing rate of the total receptor number on a cell surface per binding to the chemo-attractant receptor. So, the relation between receptors and bound receptors 13,30 is given by These two relations reduce the number of dependent variables. Finally, system (4) can be written as All the parameters involved in the system (7) are positive and summarized in the Table 2.

Non-dimensionalization.
To reduce the number of parameters, we need to do non-dimensionalization.
We choose the filter thickness L f as a length scale and one hour time T as a time scale i.e. www.nature.com/scientificreports/ cell density, chemical concentration are scaled by the initial cell density in the upper well n 0 , chemical concentration in the lower well c 0 respectively. The bound normal receptor number and bound decoy receptor are scaled by the original free normal receptor Ŵ r and decoy receptor number Ŵ d respectively. So, the dimensionless variables are The other dimensionless parameters are as follows For simplicity, we drop the tilde notation. Then, the dimensionless system can be written on the filter section i.e. on 0 < x < 1 as where K n = −D n ∂n ∂x + χn ∂ρ r ∂x − µn ∂ρ d ∂x , and on the other compartment of the Boyden's chamber i.e. on − L l < x ≤ 0 and 1 ≤ x < L u , n = n/n 0 ,c = c/c 0 ,ρ r = ρ r / Ŵ r ,ρ d = ρ d / Ŵ d . www.nature.com/scientificreports/ Initial and boundary conditions. Cell migration in response to the ligand-receptor interaction induced chemotactic sensitivity is modelled here through system (8). Hence to study cell migration in various chemotactic assays such as Boyden chamber 29 , Dunn chamber 35 , Zigmond chamber 36 etc., we need to incorporate well-defined initial and boundary conditions according to the system. Our model here is not assay specific, the boundary conditions change with different assays. In our study here, we consider the Boyden chamber's onedimensional geometry as a spatial domain of our system (8-9) which is solved using appropriate initial and boundary conditions. The initial cell number density is considered to be uniform in the upper well and zero in other compartment respectively and boundary condition follows Eq. (6). We also assume that there is no chemical in the upper well initially and the flux in both surfaces of the filters is continuous with no flux at both ends of the chamber. The initial and boundary conditions, in terms of dimensionless variables, can be summarized as followsInitial condition: Boundary condition:

Results
Sensitivity analysis. System (8) contains fourteen parameters among few parameters that were not obtained experimentally and hence we identify here significant parameters which affect the system's dynamics the most. To determine the sensitivity of these parameters, we perform a forward sensitivity analysis of a PDE model using the approach described in 37 where F (y, p, y x , y xx ) ≡ (F 1 , F 2 , F 3 , F 4 ) T are the reaction terms in the right-hand side of (8). The sensitivity index S j of the parameter p j is defined as To find the value of ∂y 1 ∂p j , we need to solve simultaneously with the system (10). More details of the method are given in the supplementary information S1. The sensitive index S j corresponding to eight important parameters namely χ, µ, β, δ, k a 2 , k d 2 , k i 2 , and Ŵ d are summarized in Table 3. The positive sign of S j implies that if the corresponding parameter p j increases then N f increases and the negative sign of S j implies that if p j increases then N f decreases for j = 1, · · · , 8. We also plot the bar diagram of S j in Fig. 2 which shows that χ, β, k d 2 and k i 2 are positively correlated to N f whereas µ, δ, k a 2 and Ŵ d are negatively correlated to N f . The most sensitive parameter is β . For the chemo-tactic parameters, the chemo-attractant parameter χ is more sensitive than the chemo-repellent parameter µ.
Numerical simulation. We perform extensive numerical simulations with the help of MATLAB software.
In our simulation, we discretize the Laplacian term and 1st order spatial derivative using central difference with mesh size h = 0.02 and time integration is performed using forward Euler scheme with t = 10 −6 . We have also checked the results with smaller choices of h and t to ensure that the obtained results are free from any numerical artifact. Results are generated upto dimensionless time t max = 5 , which corresponds to a time of 5 hours in the experimental setup or the dimensional system (7). We numerically solve system (8)(9) for the cell density n, the number of bound receptors ρ r , the number of bound decoy receptors ρ d in the filter region, and the ligand concentration c in the whole Boyden's chamber.
∂ ∂t ∂y ∂p j = ∂F (y, p, y x , y xx ) ∂p j , Table 3. Sensitive index of the parameters used in the model.  Fig. 3a. The ligand spreads throughout the whole chamber (see Fig. 3b). Chemical concentration in the upper well increases with time and at t = 5 , the ligand distribution becomes almost homogeneous. In Fig. 3c and d, the distribution of bound normal and bound decoy receptors is shown respectively. We denote distance travelled inside the filter at t = 5 by H. We determine the value of H as where n(x, t = 5) is the cell density at time t = 5. We observe that H increases with time t but when the heterogeneity of the distribution of c decreases the rate of change of H decreases too.
Variation in initial distribution. In Fig. 3, we fix n 0 = 4 × 10 5 cells/cm 3 and c 0 = 6 × 10 −8 mol/cm 3 and now we consider the initial chemical concentration c 0 and initial cell number n 0 as parameters. Since we have taken the initial distribution of cells for the dimensionless system (8) to be 1 in the upper end of the filter and 0 everywhere else in the filter, then n 0 doesn't affect the wave profile of n. But c 0 influences the system dynamics significantly. We measure distance travelled inside the filter (H) for the non-dimensional model (8)(9) for different initial ligand concentration (c 0 ) . We plot c 0 vs H in Fig. 4 and observe that initially with an increase in c 0 , H increases but after critical threshold of c T = 2.26 × 10 −8 , H decreases with an increase in c 0 .
Effect of β and δ in cell migration. Here we focus on the role of the increasing rate of normal receptors and decoy receptors β and δ in the cell migration mechanism. The distribution of the number of bound normal receptors ( ρ r ) increases with an increase in β which in turn increases the chemo-attractant sensitivity. As a result, the wave-front of the migrated cells travels much deeper inside the filter. In Fig. 5a, for β = 0.97 , the number of cells in the filter increases 85.19% more in comparison to β = 0.95. Further, an increase in δ with a fixed β , increases the distribution of decoy receptors. Therefore, cell migration is interrupted more due to additional presence of decoy receptors. For δ = 0.99 , the number of cells in the filter decreases by 19.46% in comparison to δ = 0.95 (see Fig. 5b). From Fig. 5a and b, we conclude that the chemo-attractant sensitivity of normal receptors influences the system dynamics much more in comparison to chemo-repellent sensitivity of decoy receptors.

Model validation with experimental results. Comparison of cell movement in absence of decoy receptors.
To consider the case of cell movement in absence of decoy receptors, we set all the parameters associated with decoy receptors to zero. Under this assumption, the zero initial condition of ρ d in (8d) gives only trivial solution for ρ d . This leads to the system in the filter, a three compartment coupled partial differential equations (PDE) as follows: where K n = −D n ∂n ∂x + χn ∂ρ r ∂x . System (9, 12) is subjected to same initial and boundary conditions as before. We observe that the migrated cell's wave comes deeper inside the filter in comparison to the presence of decoy receptors (see Fig. 6a). The In absence of decoy receptors, cell migration is influenced by cell random movement and the chemo-attractant sensitivity of the normal receptors. But in the case of system (8)(9), the decoy receptors inhibit the process of cell migration mechanism by its chemo-repellent sensitivity which results in a lesser number of cells in the filter. Roth et al. 8 did experiments with N9 mouse microglia cells in a 48-well micro chemotaxis chamber where  Table 2. The experimental data is provided in supplementary data. We plot the experimental data in Fig. 6b with red 'o' . In the x-axis, FL1-H scale denotes the relative concentration of CD95L. Similar pattern is also observed in our numerical simulation. We consider n 0 = 5 × 10 3 cells/cm 3 and chemical concentration c 0 as a variable. In order to match the scale with experimental results, we re-scale the initial concentration c 0 with a scaling factor 10 9 . We plot the numerical solution N f with varying c 0 in Fig. 6b with blue solid curve. Figure 6b shows that N f increases with an increase in initial ligand concentration, but after a certain threshold of ligand concentration ( c 0 [FL1-H] = 22.58 ), the relationship is reversed, with excess ligand concentration the value of N f decreases. Figure 6b shows that the numerical results are in close agreement with the experimental results in 8 .   Table 2.

Discussion
Cell migration is a directed movement of a single cell or a group of cells in response to chemo-tactic signals.
Receptor-ligand interaction is a major class of protein-protein interactions and it activates chemokine's downstream effectors which are responsible for chemo-tactic cell signalling. This chemo-tactic signalling does not always act as an attractant. A few receptors in the TNFRSF family have been identified as decoy receptors that recognize chemokines but do not induce cell migration. Roth et al. 8 did an experiment in a 48-well micro chemotaxis chamber which suggests that cell migration of gliomas cells in presence of additional decoy receptors DcR3 is significantly decreased when compared to migration in presence of only normal receptors CD95. In this work, we have extended the classic Keller-Segel model 28 to study cell migration of glioma cells in response to normal receptor CD95's chemo-attractant sensitivity and decoy receptor DCR3's chemo-repellent sensitivity. Here, we propose a mathematical model which consists of four partial differential equations that describe the spatio-temporal dynamics of cells, normal receptors, ligands and decoy receptors. We consider the Boyden chamber assay as a domain of our spatio-temporal model here but our model is not assay-specific except for the boundary conditions and we have validated this model with experimental results obtained from a 48-well micro chemotaxis chamber. We assume that initially cells and ligands are present only on the upper and lower well of the Boyden chamber respectively. The dynamics of cells, normal receptors and decoy receptors are presented upto t max = 5 in the filter section. The cell motility in our model is influenced by three factors: (1) cell's random movement, (2) chemo-attractant sensitivity of normal receptors, (3) chemo-repellent sensitivity of decoy receptor. The migrated cells move from upper well to the filter section. With the advancement of time, the spatial heterogeneity of the ligand, normal receptor and decoy receptor's distribution decreases. We observe that the initial ligand concentration plays an important role in cell migration. As seen in Fig. 4, cell migration increases with an increase in initial ligand concentration but after a certain threshold point of ligand concentration ( c T ), the relationship is reversed. One possible explanation for this phenomenon is that at high ligand concentration, the receptors on the cell surface become over-stimulated, leading to a depletion of free receptors available for binding 40,41 . This can interfere with the ability of the cells to respond to directional cues from the ligand and may even result in a loss of sensitivity to the signal. Additionally, excess ligands can potentially bind to receptors that are not involved in directing cell migration, leading to a loss of directionality 42 . This can ultimately result in a decrease in the intensity of cell movement, as the cells become less efficient at responding to the ligand signal.
Our research here successfully captures the dynamics of both the chemo-attractant activity of the normal receptors as well as the chemo-repellent activity of the decoy receptors. With an increasing rate of the normal receptors ( β ), the number of migrated cells in the filter increases whereas with an increasing rate of the decoy receptors ( δ ), the number of migrated cells in the filter decreases. We also investigate the system dynamics in the absence of decoy receptors. Figure 6a shows that the wave front curve of migrated cells is higher in the absence of decoy receptors when compared in the presence of decoy receptors. In other words, in the presence of decoy receptors, the number of migrated cells in the filter decreases. Our results support Roth's experiments 8 and we also validate our simulation results with the number of cells in the filter with changes in the initial ligand concentration. The sensitivity analysis is carried out to understand the role of important model parameters in the cell migration. Our analysis predicts that the chemo-attractant sensitivity of normal receptors is stronger than chemo-repellent sensitivity of decoy receptors.
We observe here how a cell migration model is deployed to predict the movement of glioma cells due to reaction kinetic mechanism between CD95 receptors and ligand CD95L as well as compute the distance travelled by cells inside filter. A systemic approach of exploring cell migration due to decoy receptors and studying few aspect of its movement through computational models and validating cell migration with experimental results is the strength of our study here. This work is a preliminary and a crucial step toward future in understanding the influence of decoy receptors on other mechanisms such as cell apoptosis, cell adhesion, cell signaling, and so on. Our analysis of chemotactic cell response due to receptors will have significant application and open doors to future research in quantitative understanding of immunotaxis, leukocyte migration during immune response to inflammation, and chemotactic behavior of nanoswimmers in penetration to the brain to name a few 6,7,43 .

Data availability
The datasets generated and/or analyzed during the current study are available in S2.xlsx available with the manuscript as supplementary information and also in the published article 8 .