B Cell Activation Triggered by the Formation of the Small Receptor Cluster: A Computational Study

We proposed a spatially extended model of early events of B cell receptors (BCR) activation, which is based on mutual kinase-receptor interactions that are characteristic for the immune receptors and the Src family kinases. These interactions lead to the positive feedback which, together with two nonlinearities resulting from the double phosphorylation of receptors and Michaelis-Menten dephosphorylation kinetics, are responsible for the system bistability. We demonstrated that B cell can be activated by a formation of a tiny cluster of receptors or displacement of the nucleus. The receptors and Src kinases are activated, first locally, in the locus of the receptor cluster or the region where the cytoplasm is the thinnest. Then the traveling wave of activation propagates until activity spreads over the whole cell membrane. In the models in which we assume that the kinases are free to diffuse in the cytoplasm, we found that the fraction of aggregated receptors, capable to initiate B cell activation decreases with the decreasing thickness of cytoplasm and decreasing kinase diffusion. When kinases are restricted to the cell membrane - which is the case for most of the Src family kinases - even a cluster consisting of a tiny fraction of total receptors becomes activatory. Interestingly, the system remains insensitive to the modest changes of total receptor level. The model provides a plausible mechanism of B cells activation due to the formation of small receptors clusters collocalized by binding of polyvalent antigens or arising during the immune synapse formation.


I. System bistability analysis
The bistability of the model arises due to the nonlinearities in kinases dephosphorylation (Michaelis-Menten dephosphorylation kinetics) and second order nonlinearity in receptor activation due to double phosphorylation. Each of these two nonlinearities alone, su¢ ces for the bistability. The contribution of these two nonlinearities is controlled by Michaelis constant H (the larger is the value of H the more linear is dephosphorylation kinetics) and the value of c 0 (the larger is c 0 the more linear is the receptor activation function c 0 + K 2 since always K < 1). In fact, there are additional nonlinearities in the system, that is, nonlinearity in receptor dephosphorylation, and kinase transphosphorylation. These nonlinearities are omitted for the sake of simplicity. It is natural to expect that inclusion of these nonlinearities would make the bistability region broader.
First, we con…ne ourselves to spatially uniform solutions ( = 0) and analyze the dependence of the bistability regions in the (B; H) plane on the coe¢ cient c 0 . Using Eqs.
(10),(11) (from the main text) we infer that the constant steady state solutions are given by the solutions of the system: which for c 0 = 0 simpli…es to It follows that for c 0 = 0 and for any B > 0 and H > 0 system (3) has a zero solution K = 0; R = 0. This solution is stable, because the Jacobian matrix J 0 of the vector function (R(1 K) BHK H + K ; K 2 (P R) R) calculated for K = 0; R = 0 has two negative eigenvalues equal to 1 and B. In fact, J 0 is the linearization matrix of system (10) solution" arises and from the right by line B R (H) at which the "active" solution vanishes.
Since the system behavior for c 0 = 0 is not representative, we choose c 0 = 0:01, which is approximately equivalent to the assumption that the activity of the unphosphorylated kinase is 100 times smaller than the activity of phosphorylated kinase.
Next, we set r n = 0:9, c 0 = 0:01, H = 0:1 and analyze the bistablity regions in (b; P ) plane for three values of . As shown in Figure S2 the bistability range in (b; P ) plane for = 3 is almost identical to the bistability range for the in…nite di¤usion, = 0. This shows that the analysis presented in Figure S1 can serve as a rough reference also for the …nite but large di¤usion 2 .

II. Conditions for activatory traveling wave propagation
The traveling waves observed numerically both in the cytosolic and the membrane models are not the usual plane traveling waves. In the cytosolic model, at a given moment of time, the values of K and R are not uniform along the radius direction. In both models the curvature of the wave front is very important -especially when the activation starts from receptors cluster occupying a small fraction of the cell membrane. Here, basing on the membrane model, we will discuss the necessary conditions for activatory traveling wave propagation, then we will provide some numerical examples ( Figures S5 and S6) for the cytoplasmic model for which the analysis is more complicated.
In the membrane model the equation for the active kinase concentration K may be written as The …rst term of the right-hand-side of Eq. (4) is associated with the wave front curvature.
Presence of the curvature term causes that the wave front propagation velocity depends on . As shown in Figure 7 (in the main document) in the main document, the velocity of activatory front increases with . Without the curvature term the "membrane model" (for a = 1; P = 1) reads and may be considered on the in…nite domain z 2 ( 1; 1). It has thus traveling wave solutions K (t; z) = K(z ct); R (t; z) = R(z ct): In this case, the propagation velocity c is c = 1 f (b; H; c 0 ). For …xed H = 0:1; c 0 = 0:01, the system is bistable for b 2 [6:73; 27:1] and the bistability range is independent of . The critical value of b (b crit ) corresponds to standing wave solutions with c = 0. For b < b crit the activatory traveling waves may propagate, while for b > b crit , traveling waves propagate in the opposite direction -thus they can be considered de-activatory. For the reduced system (5)-(6) b crit = 20:5.
The sign of the curvature term 2 cot( )@K=@ , for monotone traveling wave pro…le (@K=@ < 0, see Figure 7 in the main document), is determined by the sign of the function cot( ). For < =2 presence of the curvature term decreases @K=@t and thus traveling wave velocity (assuming that the wave propagates from = 0 to = ), while for > =2 the curvature term increases the propagation velocity. For small values of , the contribution of the curvature term may dominate the right-hand-side of Eq. (4) and change the sign of the wave front propagation, or prevents the front from leaving the receptors cluster area in which R is large. One may thus expect that the critical value of b, below which the activatory wave will leave the receptor cluster will be smaller than the critical value of b found for plane front, b crit = 20:5. In order to demonstrate the curvature e¤ect we consider the membrane model with the following receptor distribution P ( ) = 1 + i + 1 2 i (1 + cos ) i F and calculate b crit for two values of i, i = 1000 and i = 10000 and three values of . The large value of F = 0:1 implies the local system activation (see Figure S4) but does not imply that activatory wave front will propagate over the area where the receptors density is close to 1. For i = 1000 we obtained the critical values of the parameter b, as 16:3, 17:6 and 19 for = 1, = 3 and = 10 respectively. For i = 10000, b crit is equal 14:4, 15:8 and 17:8 for = 1, = 3 and = 10 respectively. The curvature e¤ect, i.e. the di¤erence between the critical value of b, b crit = 20:5, obtained for plane traveling waves and b crit for the membrane model (5)-(6), increases with increasing i (i.e. initial front curvature) and with decreasing . The curvature e¤ect causes that the clusters of very small area are not activatory.