Trait-space patterning and the role of feedback in antigen-immunity coevolution

Coevolutionary arms races form between interacting populations that constitute each other’s environment and respond to mutual changes. This inherently far-from-equilibrium process finds striking manifestations in the adaptive immune system, where highly variable antigens and a finite repertoire of immune receptors coevolve on comparable timescales. This unique challenge to the immune system motivates general questions: How do ecological and evolutionary processes interplay to shape diversity? What determine the endurance and fate of coevolution? Here, we take the perspective of responsive environments and develop a phenotypic model of coevolution between receptors and antigens that both exhibit cross-reactivity (one-to-many responses). The theory predicts that the extent of asymmetry in cross-reactivity is a key determinant of repertoire composition: small asymmetry supports persistent large diversity, whereas strong asymmetry yields long-lived transients of quasispecies in both populations. The latter represents a new type of Turing mechanism. More surprisingly, patterning in the trait space feeds back on population dynamics: spatial resonance between the Turing modes breaks the dynamic balance, leading to antigen extinction or unrestrained growth. Model predictions can be tested via combined genomic and phenotypic measurements. Our work identifies cross-reactivity as an important regulator of diversity and coevolutionary outcome, and reveals the remarkable effect of ecological feedback in pattern-forming systems, which drives evolution toward non-steady states different than the Red Queen persistent cycles.


I. INTRODUCTION
Highly variable antigenic challengers, such as fast evolving viruses and cancer cells, are rapid in replication and abound with genetic or phenotypic innovations [1,2], thus managing to evade immune recognition. On the reciprocal side, the host immune system adjusts on the fly the clonal composition of its finite receptor repertoire to recognize the altered versions of the antigen. As a result, coevolutionary arms races between antigen and immunity may endure through an individual's lifetime [3].
In this Red Queen [4] scenario, antigen and receptor populations constitute each other's responsive environment and are mutually driven out of equilibrium: specific immune receptors prey on matching antigens and hence alter both the composition and overall abundance of antigens, which in turn modifies selective pressures on distinct receptors, thus causing reorganization of the repertoire, and vice versa. Consequently, neither population has enough time to equilibrate and yet they mutually engage in a dynamic balance. In this sense, the Red Queen state represents a nonequilibrium steady state [5,6]. Then, the question is whether alternative evolutionary outcomes characteristic of nonsteady states occupy a larger * shenshen@physics.ucla.edu volume of the state space of coevolving systems than does the Red Queen state.
Recent progress has been made toward understanding various aspects of coevolutionary dynamics in antigen-immunity systems [7][8][9][10][11][12][13][14][15], ranging from antibody evolution against HIV and influenza viruses to evolution of tumors and bacterial phage under host immunity. Yet, we are still short of insights into certain fundamental questions: How do receptor repertoire and antigen ensemble mutually organize, when ecological and evolutionary dynamics occur on comparable timescales? What governs the persistence and outcome of mutual adaptation?
In existing generic models where both ecological and evolutionary processes are considered, a separation of timescales is often assumed so that the fast dynamics is slaved to the slow one (reviewed in [16]). In cases where timescales are not treated as separated [17][18][19][20][21][22][23][24], the feedback between changes in diversity and population dynamics tends to be ignored. The goal of this paper is to consider inseparable timescales and at the same time account for feedback effect in order to address the questions raised above.
Specifically, we develop a phenotypic model, based on predator-prey interactions between coevolving immune receptors and antigens, that combines evolutionary diversification and population dynamics. By formulating an ecological model in a trait space, we describe coevolutionary changes in the distribution of trait values and trait-dependent predation in the same framework. Importantly, this allows us to study the stability of speciation (pattern formation in the trait space) and its impact on the persistence of coevolution. Our model abstracts the key features of adaptive immunity: antigens and receptors move (due to trait-altering mutations) and behave like activators and inhibitors that react through predation; both antigens and receptors are cross reactive (one receptor recognizes many distinct antigens and one antigen is recognized by multiple receptors), this flexibility in recognition stems from structural conservation of part of the receptor/antigen binding surface [25] and provides an enormous functional degeneracy [26] among distinct immune repertoires; there is no preexisting fitness landscape for either population so that selection pressures are owing purely to predation.
The theory predicts, counterintuitively, that simultaneous patterning in coevolving populations can emerge solely from asymmetric range of activation and inhibition in predator-prey dynamics, without a need for severely large differences in their rate of evolution [27] (aka mobility in their common phenotypic space), thus representing a Turing mechanism distinct from the classic one. This surprising result can be understood from an intuitive picture: colocalized clusters of antigens and receptors form in the trait space when the "inhibition radii" of adjacent receptor clusters overlap so that inhibition of antigen is strongest in-between them; whereas alternate clusters emerge if the "activation radii" of neighboring antigen clusters intersect because then activation of receptor is most intense in the midway. Biologically, receptor activation and antigen inhibition are distinct processes: the assumed asymmetry in reaction range reflects potential distinction between the ability of antigens to induce protective immune responses (immunogenicity) and the ability to interact with the product, such as antibodies, formed by a response (antigenicity). In fact, this discrepancy between antigenicity and immunogenicity has been known for long [28] and demonstrated for both natural and synthetic antigens [29][30][31].
We show that as asymmetry in cross reactivity varies, transitions between qualitatively distinct regimes of ecoevolutionary dynamics seen in nature would follow, including persistent coexistence, antigen elimination, and unrestrained growth. While competitive interactions, whether direct or mediated by resource competition, are known to elicit patterns in a population [32][33][34][35], an interesting outcome of our analysis is that mutual feedback between dynamic patterns of antigens and receptors can drive the arms race off balance. Given sufficient asymmetry, spontaneous oscillations in Turing patterns precede antigen extinction, whereas uncontrolled antigen growth follows the formation of alternate quasispecies, as ineffective receptors exhaust the limited immune resources; these measurable features may serve as precusors of the offbalance fates.
Many theoretical studies have considered adaptation to time-varying environments with prescribed environmental statistics [36][37][38][39][40][41][42]. This work makes a step toward a theory of coevolution from the perspective of responsively changing environments (mutual niche construction [43] in ecological terms), highlighting the role of feedback in driving evolution toward novel organization regimes and nonsteady states. As new genomic and phenotypic methods are developed to better characterize antigenic [44][45][46] and immunological [47][48][49][50] landscapes as well as bidirectional cross reactivity [51], the predictions for repertoire composition and coevolutionary outcome derived from this study can be compared with high-throughput profiling of coevolving immune repertoire and antigen ensemble in humans [52][53][54].

II. MODEL
A finite repertoire of immune receptors that collectively cover the antigenic space while leaving self-types intact is conceivable, if the distribution of potential threats is predetermined [55][56][57][58]. Given a fixed distribution of pathogenic challenges, competitive exclusion is shown to drive clustering of cross-reactive receptors [57]. In coevolution, however, antigen distribution is no longer preset but responds to reorganization of receptors. In addition, cross reactivity is bidirectional: not only can a receptor be activated by a range of distinct antigens, but an antigen can be removed by a variety of receptors. Then, can predation lead to simultaneous clustering of antigens and receptors in their common trait space? If so, are such patterns stable? Would the concurring patterns interact to affect population dynamics? To answer these questions, we consider a dynamical system of activators and inhibitors representing antigens and receptors, which diffuse in a shared phenotypic space and react through predator-prey interactions. Population densities of antigens A(⃗ x, t ) and receptors B(⃗ x, t ) evolve according to Here, D 1 and D 2 denote isotropic diffusion constants of antigens and receptors, respectively, that mimic the rates of traitaltering mutations. Other forms of jump kernels do not change qualitative results, see Supplemental Material Ref. [66]. Antigens self-replicate at rate λ 1 whereas receptors spontaneously decay at rate λ 2 . Receptors inhibit antigens with an intrinsic rate α 1 and grow at rate α 2 upon activation; R inh denotes the range of receptors that can inhibit a given antigen, while R act represents the range of antigens by which a receptor can be activated (Fig. 1). In real systems, there is likely a distribution of reaction range; we assume a single value to simplify analysis. The term B in corresponds to a small influx of lymphocytes that constantly output from the bone marrow and supply nascent receptors; without stimulation, receptors are uniformly distributed at a resting concentration given by B in /λ 2 . We choose the lifetime of receptors λ −1 2 as the time unit and the linear dimension L of the phenotypic space as the length unit. To account for the discreteness of replicating entities and hence avoid unrealistic revival from vanishingly small population densities, we impose an extinction threshold; antigen or receptor types whose population falls below this threshold are considered extinct and removed from the system.
In the spirit of Perelson and Oster [59], we think of receptors and antigens as points in a high-dimensional phenotypic space, whose coordinates are associated with physical and biochemical properties that affect binding affinity. We assume that the strength of cross-reactive interaction only depends Schematic of antigen-receptor interaction with asymmetric range of inhibition and activation in the phenotypic space. (a) R inh > R act : the receptor (blue Y shape) is not activated by the antigen (red flower shape) but nevertheless inhibits it. (b) R act > R inh : the antigen activates the receptor but is not subject to its inhibition. Lower row: in addition to predation (black arrows; blunt for inhibition, acute for activation), antigens self-replicate (red arrow) whereas receptor-expressing cells spontaneously decay (blue arrow pointing to an empty set symbol) in the absence of stimulation. If a finite carrying capacity of receptors θ 2 is explicitly considered, self-inhibition will also be present [Figs. 4(c), 6(b), and 6(c)].
on the relative location ⃗ r = ⃗ x − ⃗ y, of receptor and antigen in this space, as characterized by the nonlocal interaction kernels S 1 (|⃗ r|; R inh ) and S 2 (|⃗ r|; R act ). Close proximity indicates a good match between the binding pair leading to strong interaction, whereas large separation translates into weak affinity and poor recognition.
Importantly, cross reactivity is not necessarily symmetric as typically assumed; differences in biophysical conditions among other factors may well render disparate criteria for antigen removal and receptor activation [60,61], i.e., R inh ̸ = R act . For instance, removing an antigen may only require modest on rate (wide reaction range, large R inh ) of multiple receptors that together coat its surface, which boils down to multivalent binding and cross linking at thermal equilibrium. Whereas activating an immune cell expressing a unique type of receptors can demand lasting antigen stimulation hence small off rate (close match of shape, small R act ), so that downstream events leading to a response can finish. How this asymmetry impacts coevolution is our focus.

A. Phases under local predator-prey interactions
This reaction-diffusion system [Eq. (1)] presents a homogeneous fixed point of population densities A s = λ 2 /(α 2 2 ) and B s = λ 1 /(α 1 1 ), where 1 = d⃗ r S 1 (|⃗ r|; R inh ) and 2 = d⃗ r S 2 (|⃗ r|; R act ) are, respectively, the shape-space volume of the "inhibition sphere" centered at an antigen and that of the "activation sphere" surrounding a receptor. When receptor-antigen interactions are local, depending on the ratio of the rates λ 1 /λ 2 and diffusivity D 1 /D 2 , coevolving populations exhibit two main phases within the chosen parameter range [ Fig. 2(a)]: antigen early extinction (colored region) and persistence (white region); the latter divides into two subphases, steady traveling waves (upper) and uniform  Phases in a 1D reaction-diffusion system under local predator-prey interactions. (a) Phase diagram on the plane spanned by the ratio between diffusion constants D 1 /D 2 and that between birth and death rates λ 1 /λ 2 of antigens (activators) and receptors (inhibitors). Dynamics start from a local dose of antigens and uniform receptors. The early extinction phase is color coded for the logarithm of the inverse time to antigen extinction. The persistence phase (blank) divides into a propagating wave state (upper) and a uniform coexistence state (lower). Insets show typical kymographs in each subphase, red for antigen and blue for receptor; the upper pair corresponds to the filled circle at λ 1 /λ 2 = 200, D 1 /D 2 = 10 −2 , and the lower one corresponds to the open circle at λ 1 /λ 2 = 10, Corresponding phase plots are shown on the right; vertical dashed lines indicate the extinction threshold. B in = 10, α 1 = 10 −3 , α 2 = 10 −4 . coexistence (lower); as seen in typical kymographs of the one-dimensional (1D) density fields (insets), starting from localized antigens and uniform receptors.
Extinction is expected when antigens replicate fast (large λ 1 /λ 2 ) but mutate slowly (small D 1 /D 2 ): after a brief delay during which antigen reaches a sufficient prevalence to trigger receptor proliferation, receptors rapidly expand in number and mutate to neighboring types; the pioneer receptors stay ahead of mutating antigens and eliminate them before escape 033164-3 mutants arise. Once antigen is cleared, the receptor population regresses to the resting level [ Fig. 2(b), upper panel]. With sufficiently high replication rates, faster mutation allows antigen mutants to lead the arms races against receptors resulting in a persistent evolving state, a traveling wave Red Queen state, similar to that shown in a recent model of influenza evolution under cross immunity [62]. Interestingly, as λ 1 /λ 2 increases, while the rate of extinction (color coded in the phase region) increases, a smaller D 1 /D 2 is needed for transition to the traveling wave state. On the other hand, at modest λ 1 /λ 2 , a uniform coexistence phase is reached following population cycles dampened by mutation [ Fig. 2(b), lower panel]. Under local interactions, this homogeneous fixed point is stable to perturbation and does not support spontaneous antigen speciation (i.e., breakup of a continuum into fragments in the shape space). Thus, in what follows, we start from this uniform steady state and introduce the key ingredient, asymmetric nonlocal interaction, to show how it drives spontaneous organization.

B. Simultaneous patterning under asymmetric cross reactivity
The analogy between antigen-immunity interaction and predation has been made before [8,15,63,64]; however, spontaneous speciation has not been described yet. On the other hand, for general activator-inhibitor systems, in the physical space, Turing patterns can emerge, either from demographic stochasticity which prevents the system from reaching its homogeneous fixed point [35,65] or, more classically, from prohibitively large differences in diffusivity between the autocatalytic and the inhibitory reactants [27]. Here, using a simple phenomenological model accounting for cross reactivity [Eq. (1)], we show that coevolutionary speciation is possible, without requiring any of the aforementioned patterning mechanisms.
To identify the onset of patterning instability, we perturb the uniform stationary state (A s , B s ) with nonuniform varia- where ⃗ k is the wave vector of a spatial Fourier mode. The Turing instability occurs when the most unstable mode, i.e., the critical mode with a wave number k c , begins to grow while all other modes decay. This condition, Re[ω(k c )] 0, corresponds to the following inequality in the n-dimensional recognition space: whereŜ 1 (k) andŜ 2 (k) are the Fourier transform of the interaction kernels. It immediately follows that Turing instability in our system is purely driven by asymmetric nonlocal interactions and independent of diffusion: if the kernels were symmetric, i.e., S 1 (r; R inh ) = S 2 (r; R act ), the righthand side of Eq. (2) can never be positive and hence patterns do not develop; on the other hand, when D 1 D 2 = 0, the patterning condition is most readily satisfied, implying that diffusion is not necessary. In fact, the commonly assumed Gaussian kernel represents a marginal case which does not robustly warrant instability [33]. Instead,Ŝ(k) < 0 is guaranteed if the strength of interaction decreases steeply with increasing separation across the edge of the interaction range, see Supplemental Material Ref. [66]. For sim-plicity, we assume step-function kernels S 1 (r) = (R inh − r) and S 2 (r) = (R act − r). Accordingly, under a modest extent of asymmetry, i.e., γ ≡ (R inh − R act )/(R inh + R act ) ≪ 1, the pattern-forming condition can be explicitly expressed in terms of γ , see Supplemental Material Ref. [66]: or, equivalently, Here, C n is a constant that only depends on the dimension n of the shape space. Rapid increase of C n with n ( Fig. S3 in Supplemental Material Ref. [66]) indicates that stable uniform coexistence extends to stronger asymmetry as the phenotypic space involves higher dimensions. Therefore, under sufficient asymmetry, a continuum of antigen (receptor) types spontaneously segregates into species-rich and speciespoor domains with densities on either side of A s (B s ). The spacing between adjacent antigen or receptor density peaks, i.e., the pattern wavelength λ ≃2π /k c , is modestly larger than the sum of activation and inhibition radii [ Fig. 3(a)] due mainly to asymmetry and slightly to diffusion. Note that the minimum level of asymmetry required for patterning decreases with increasing range of cross reactivity as γ c ∼ (R inh + R act ) −2 [Eq. . This seemingly counterintuitive behavior can be explained by a rather general mechanism. When 2R inh > λ, the "inhibition sphere" of an antigen may enclose adjacent receptor density peaks. As a result, locations between the peaks, where the actual receptor density B(x) (blue solid line) is in fact the lowest, are instead the worst positions for antigens to be in because the effective receptor density field acting on antigens at position x, B eff (x) = x+Rinh x−Rinh B(y)dy (blue dashed line), is maximal when x is right amid receptor peaks. Thus, the antigen distribution winds up tracking the receptor distribution (yellow bar, left panel). Conversely, when 2R act > λ, the "activation sphere" of a receptor may encompass adjacent antigen peaks; the stimulation for receptor replication is strongest in-between the peaks, according to the effective antigen densities A eff (x) = x+Ract x−Ract A(y)dy (red dashed line). Consequently, receptors view antigens as most concentrated in positions where they are actually least prevalent (red solid line). Therefore, depending on whether R inh or R act is larger, colocalized or alternate distributions result, which reflect a mismatch between the actual distribution and the effective one seen by the apposing population. In what follows, we show that distinct spatial phase relations between mutual distributions will lead to drastically different pattern dynamics and evolutionary outcomes. Intriguingly, concentrations and patterns evolve via three distinct stages: uniform steady state, stationary pattern, and oscillatory pattern. Right after copatterns spontaneously emerge from the homogeneous steady state, concentration changes in both populations are observed: in the absence of homeostatic constraints [Figs. 4(a) and 4(b)], antigen abundance (red) shifts downward whereas receptor prevalence (blue) shifts upward. This is unanticipated because patterning instability in a density field is not expected to alter the overall abundance: growing unstable modes merely redistribute den-sities in space without changing the average concentration. This appears to break down when patterns develop in two interacting density fields. In fact, the most unstable modes (with wave number k c ) from both populations couple and modify the zero modes, resulting in a shift in mean population densities.
Colocalized and alternate quasispecies. A weakly nonlinear analysis close to the critical point quantitatively captures both the phase relation between patterns and the shift in overall abundances (Fig. 5). For analytical tractability, we perform the calculation in 1D, see Supplemental Material Ref. [66]. Below, we only stress the essential results.
Close to the patterning transition 1 is the critical diffusion constant of antigen and ϵ is small and positive, we seek stationary solutions of the form where the deviation w = (u(x), v(x)) T from the homogeneous steady state (A s , B s ) T is expanded in powers of ϵ 1/2 to the second order: The saturated amplitude A of the perturbation is determined by the amplitude equation at the order of ϵ 3/2 . The spatial phase difference between the leading pattern modes in antigen and receptor populations, u (1) 1 and v (1) 1 , respectively, can be found from which implies that ξ ∝ sin ( 2π λ (Ract+Rinh ) 2 (1 − γ )) = sin (π 1−γ λ/(Ract+Rinh ) ), with λ being the pattern wavelength and γ = (R inh − R act )/(R inh + R act ). It immediately follows that when γ < 0 (i.e., R act > R inh ), 1 − γ > λ/(R act + R inh ) > 1 [see Fig. 3(a)], thus ξ < 0; when γ > 0 (i.e., R inh > R act ), 1 − γ < 1 < λ/(R act + R inh ), thus ξ > 0. Therefore, the spatial patterns of antigen and receptor distributions are either in phase (ξ > 0) or out of phase (ξ < 0), purely determined by the sign of asymmetry [ Fig. 5(a)]. This provides rigor to the intuitive argument we made earlier in relation to Fig. 3(c).

indicating that shift in abundance indeed results from coupling between simultaneous Turing modes.
Dynamic transients. A further surprise comes at longer times: the stationary copatterns are only metastable. Soon after abundance shift takes place, instability starts to grow, visible as increasingly strong oscillations that eventually drive the antigen population to pass below the extinction threshold [ Fig. 4(a), top panel]. By perturbing around the abundanceshifted stationary patterns, we indeed identify a growing oscillatory instability from the dispersion relation of linearized dynamics, see Supplemental Material Ref. [66]. The interrupted oscillation amplitudes at later times arise from asynchronous extinction of local antigen clusters [ Fig. 4(a), middle panel]. Note that this late extinction phase only occurs to colocalized population densities, i.e., when R inh > R act .
Upon interchange of R inh and R act [ Fig. 4(b)], pattern evolution exhibits different features: as some antigen clusters go extinct as a result of oscillatory instability, neighboring clusters migrate to these just vacated sites, where receptors decay due to a lack of stimulation and delay in response, thus locally and temporarily evading immune inhibition. These surviving clusters then go through successive branching events (i.e., widening then splitting), forming a treelike structure over time. The coevolving receptor population drives the branching and subsequently traces the newly formed branches (see movie a in Supplemental Material Ref. [66] for pattern dynamics). Such coevolutionary speciation, enabled by mutation, persists for extended periods of time so that it effectively overcomes the growing oscillations and maintains antigen at modest prevalence indefinitely. Note that the persistent ramifying pattern only emerges from alternate density peaks, i.e., when R act > R inh .
Finite repertoire. Is there a chance that antigen population can achieve a global escape from immune control, as is often envisaged as a catastrophic failure? This does happen as soon as we turn on a sufficiently strong homeostatic constraint on receptor abundance [ Fig. 4(c)]. Considering global homeostasis, the decay rate of receptors becomes λ 2 B(x, t )(1 + 1 0 B(y, t )dy/θ 2 ), which now includes an additional contribution from the global constraint characterized by carrying capacity θ 2 (see Supplemental Material [66], Sec. IV, for discussion). Interestingly, reducing the immune capacity appears to alter the nature of the instability [ Fig. 6(c), Fig. S9 in Supplemental Material Ref. [66]]: a critical value of θ 2 marks the transition from supercritical bifurcation (yellow region), where nonlinearity acts to saturate the growth of the perturbation, to subcritical bifurcation (red region) where higher-order processes have to intervene for stabilization. The latter corresponds to the antigen escape phase [ Fig. 4(c)]: an unrestrained growth indicates a loss of immune control.
Higher dimension. Similar progression of patterns and population dynamics in distinct regimes is also seen in two dimensions (2D) starting from the uniform steady state (movies b to d in Supplemental Material Ref. [66]). An analogous "branching" scenario in the persistence phase is particularly intriguing: antigen droplets deform and migrate to neighboring vacant loci and resist elimination. Oscillations of dense spots in both populations resemble the "twinkling eyes" pattern proposed for synthetic materials. It has been suggested [67] that oscillatory patterns can arise in a system consisting of two coupled reaction-diffusion layers, one capable of producing Turing patterns while the other supporting Hopf instability. Distinct from these built-in mechanisms, instabilities in our system are self-generated: interacting populations spontaneously fragment in the trait space, and the resulting Turing modes resonate in space, leading to abundance shift and subsequent growing oscillations.
Phase diagram. To stress the role of asymmetric cross reactivity in governing the diverse behaviors, we present phase diagrams on the (R act , R inh ) plane (Fig. 6). Without homeostatic constraints [θ 2 = ∞, Fig. 6(a)], patterns form above the critical asymmetry marked by solid lines that are symmetric about the diagonal; the enclosed patternless phase (light blue region) corresponds to stable homogeneous coexistence like for local interactions. On the R inh > R act side, the late extinction phase (blue) transitions to the persistent patterned phase (yellow) at a boundary (dashed line) determined by tracking the prevalence trajectories until t = 100; longer tracking time would expand the extinction phase.
Under a finite carrying capacity [θ 2 = 3 × 10 5 , Fig. 6(b)], two major changes occur: First, the pattern-forming region is no longer symmetric but expands on the R act > R inh side toward the diagonal. Second, the antigen escape phase (red region) emerges at the small-R inh large-R act corner; the phase boundary corresponds to the transition between supercritical and subcritical bifurcations in the amplitude equation. The escape phase enlarges as the carrying capacity diminishes. Thus, our model predicts expansion of the antigen escape phase with age, owing to diminishing counts of renewable lymphocytes [68]. The regime of persistence with pattern (yellow), irrespective of the homeostatic constraint, differs between the flanks; while oscillations occur on both sides off the diagonal, antigen branching [ Fig. 4(b)] only appears when R act > R inh , manifesting the potential for evasion.

IV. DISCUSSION
Environment becomes a relative concept in light of coevolution. We present a general model of mutual organization between continuous distributions of antigens and receptors that interact cross reactively. In a shared phenotypic space, the receptor repertoire and antigen population constitute each other's environment and adapt to mutually constructed fitness seascapes. This phenomenological approach allows us to describe the interplay between ecological and evolutionary processes that do not separate in timescales, thus revealing a variety of dynamic transients observed in nature, such as antigen extinction, chronic persistence, and unrestrained growth until saturation.
We propose that the transient nature of host-pathogen coevolution could, at least in part, stem from distinct conditions for receptor activation and antigen inhibition. On the one hand, the ability of antigens to be recognized by the immune system, i.e., antigenicity, can be reduced to the level of chemistry and measured by in vitro lymphocyte proliferation and cytokine production. On the other hand, the ability to induce protective immunity, i.e., immunogenicity, depends on complex interactions with various elements in the host immune system, thus demanding immunization studies in vivo. Indeed, experiments have demonstrated for diverse pathogens that strong antigenicity does not guarantee protection and vice versa; this lack of correlation has posed significant challenges to vaccine design [29][30][31].
Our simple model accounts for this intrinsic asymmetry and predicts its influence on antigen-immunity coevolution. While it might be intuitive that under reciprocal cross reactivity, antigen and receptor populations simultaneously fragment in the phenotype space (Fig. 3), more surprises come after the copattern emerges (Fig. 4): When two distributions are in phase (R act < R inh ), spatial resonance between the lowest Turing modes precedes growing oscillations in the overall abundance, driving antigens to extinction; when apposing populations are out of phase (R act > R inh ), strong homeostatic constraints on immune cells alter the nature of pattern instability from supercritical to subcritical, leading to uncontrolled growth. The intuitive picture is, when R act < R inh , antigens are inhibited by receptors that they do not activate and hence fail to evade immune attack; when R act > R inh , receptors are activated by antigens that they cannot inhibit, thus, under resource limits, an increasingly weaker defense results. Such multistage patterning and its feedback to population dynamics, triggered by asymmetric nonlocal interactions, is a qualitatively unique phenomenon, clearly different from speciation due to competitive exclusion in a single population. Our predictions are supported by experiments: strong oscillations in antigen abundance prior to crash to extinction have been seen in viral evolution within humans and attributed to cross-reactive antibody response [69], whereas strategies of distracting immune attention are indeed used by many viruses that create a vast excess of defective particles than functional ones [70].
These predictions can potentially be tested by tracking both the pathogen load and diversity history via high-throughput longitudinal sequencing of receptors and antigens [3,[52][53][54]. In addition, phenotypic assays for binding and neutraliza-tion [51] can inform the extent of asymmetry. Combining these two sets of experiments in different individuals would allow to correlate the degree of asymmetry with evolutionary outcomes.
Our results also suggest that the immune system may have evolved to exploit the asymmetry between activation and inhibition by differentiating these processes physically and biochemically. A remarkable example is affinity maturation of B lymphocytes [71] in which rapid Darwinian evolution acts to select for high-affinity clones: immature B cells are trained in lymphoid tissues where antigens are presented in a membrane form and decline in availability; fierce competition for limited stimuli thus provides a sustained selection pressure that constantly raises the activation threshold, i.e., decreasing R act . In contrast, mature B cells then released into circulation encounter soluble antigens at higher abundance, corresponding to R inh > R act . As a result, enhanced asymmetry between conditions for immune stimulation and antigen removal facilitates elimination of pathogens. Conversely, pathogens evolve immunodominance [72] and make fitness-restoring mutations [73] that increase R act and decrease R inh , both of which aid in evasion.
Another implication of this study is persistent coevolution, often pictured as an asymptotic state, can only be sustained when asymmetry is not too strong. It might be favorable if asymmetry stays near the edge between persistence and imbalance, which adjusts to the tension between the need for defense against foreign pathogens (γ > γ c ) and that for tolerance toward benign self-tissues (|γ | γ c ). Interestingly, critical asymmetry γ c increases with the number of phenotypic dimensions (Fig. S3 in Supplemental Material Ref. [66]), suggesting that dynamic balance could be easier to maintain for more complex antigens. Our analysis also reveals other stabilizing factors for coexistence, see Supplemental Material Ref. [66], including a stronger influx of naive cells, larger jump sizes in trait values, and longer-tailed interaction kernels.
Because the present model of antigen-immunity coevolution is a sufficiently abstract one, having properties which seem quite robust and independent of the details of predation, we expect that the results and predictions are relevant for a wide range of coevolving systems, including cancer cells and T lymphocytes, embryonic tissues and self-reactive immune cells, as well as bacteria and bacteriophage. This model can be adapted to be more biologically faithful, e.g., by incorporating preexisting antigenic landscapes, taking rates to be age dependent, and treating cross reactivity as an evolvable character.
Stochasticity arising from demographic noise does not change qualitative model behaviors in all regimes, see Supplemental Material Ref. [66]. Albeit not required for pattern formation, stochastic fluctuations appear to speed up instability growth thus accelerating antigen extinction (Fig.  S12 in Supplemental Material Ref. [66]); this observation and other effects of demographic noise will receive a careful analysis in future work.
We hope that this work proves useful in providing a framework for understanding and testing how cross-reactive interactions, ubiquitous and crucial for biological sensory systems,  can lead, in part, to the generation, maintenance, and turnover of diversity in coevolving systems. More broadly, our work provides the basis for a theory of evolution in responsively changing environments, highlighting that ecological feedback in pattern-forming systems can yield dynamic transients and drive evolution toward nonsteady states that differ from the Red Queen persistent cycles.

ACKNOWLEDGMENTS
We thank Sidney Redner and Paul Bressloff for enlightening discussions. S.W. gratefully acknowledges funding from the Dean of Physical Sciences at University of California, Los Angeles.

APPENDIX
Below, we provide additional description of the main steps in our analysis. The three subsections are ordered in accordance with the successive instabilities that give rise to the multistage progression of trait-space patterns and associated population dynamics.

Linear instability of the uniform steady state
To identify the onset of patterning instability, we linearize the equation of motion [Eq. (1) in the main text] around the homogeneous fixed point (A s , B s ) and work in the Fourier space. Defining A(x, t ) whereŜ 1 (k) andŜ 2 (k) are the Fourier transforms of interaction kernels S 1 (r) and S 2 (r), respectively. Solving this characteristic equation gives the dispersion relation, i.e., the linear growth rate of the Fourier modes: Turing instability occurs when the least stable mode (with a wave vector k c ) begins to grow, namely, where the wave number of the critical mode can be determined by ∂ k ω| k=kc = 0. This gives the pattern-forming condition [Eq.

Weakly nonlinear stability analysis: Amplitude of patterns
In this section we lay out the procedure of amplitude expansion [74,75] that yields an approximate description of the patterns formed.
Close to the transition, we introduce a small positive parameter ϵ = (D ⋆ 1 − D 1 )/D ⋆ 1 . Thus, terms in the equation of motion can be collected into a linear part evaluated at the transition L ⋆ w, and the rest g ϵ (w), for which the evolution operator explicitly depends on ϵ: Here, w = (u(x), v(x)) T denotes the deviation from the homogeneous steady state.
Note w n ∝ A n , where A represents pattern amplitude; the closer to the onset of instability, the smaller the saturated amplitude of the perturbation. Substituting this expansion into the equation of motion yields order by order (explicit expressions can be found in Supplemental Material Ref. [66]) L ⋆ w w w 2 = F F F (w w w 1 ), O(ϵ 3/2 ), L ⋆ w w w 3 = G G G(w w w 1 , w w w 2 ). The first equation effectively recovers the linear theory. For nontrivial solutions to these equations to exist, certain conditions need to be met. Such solvability condition associated with the third equation serves to determine the pattern amplitude A. It reads as where ρ † is the nontrivial kernel of the adjoint operator L ⋆ † , satisfying L ⋆ † ρ † = 0. This orthogonal condition results in the stationary amplitude equation (or Stuart-Landau equation) a 0 A − a 1 A 3 = 0, where a 0 is the linear growth rate and a 1 the Landau parameter, both of which depend on system parameters. The formation of steady patterns with a finite amplitude requires a 0 × a 1 > 0, corresponding to supercritical bifurcation where nonlinearity acts to saturate the growth of linearly unstable modes. In contrast, a 0 × a 1 < 0 leads to subcritical bifurcation, indicating that even higher-order nonlinearities are required to stabilize the pattern.
This method also applies when homeostasis is considered. Transition from supercritial to subcritical bifurcation occurs when the carrying capacity of immune receptors (predators) falls below a critical value [ Fig. 6(c) and Fig. S9 in Supplemental Material Ref. [66]).
As such, we obtain an approximate solution of the patterned state close to transition [Eq. (5) in the main text], which shows nice agreement with the numerical solution of the equation of motion (Fig. 5).

Pattern stability
To determine the stability of the patterned state, we introduce a small perturbation δw to the steady pattern w s obtained in the last section, i.e., w = w s + δw, so that the perturbation evolves according to ∂ t δw = L s δw, where L s stands for the linearized evolution operator evaluated at w s . Again, we expand the perturbation in Fourier 033164-9 modes δw = ∞ −∞ δw q e iqx dq. Since the operation of L s contains modes k c and 2k c [according to Eq. (5)], it couples q mode with those having wave numbers q ± k c , q ± 2k c , etc. Consequently, unlike perturbations to the homogeneous state that lead to decoupled characteristic equations for individual modes, now modes that differ by multiples of k c are all coupled and evolve together according to a matrix equation In principle, the state vector would be infinite dimensional, containing modes of perturbations associated with wave numbers q, q ± k c , q ± 2k c , . . . . Correspondingly, the evolution matrix M would also be infinite, forbidden from being solved. Thus, truncation of the cascade is needed to make progress.
After truncating the matrix to a modest rank, we can proceed by numerically solving for the eigenvalues of the characteristic matrix. Then, like in usual stability analysis, the leading eigenvalue m 1 (q) describes the long-term growth of the perturbation. Close to the transition, m 1 (q) solved as such agrees well with the linear growth rate of perturbations observed in numerical solutions of the original equation of motion (Fig. S7 in Supplemental Material Ref. [66]). In particular, this analysis captures the growing oscillations in the patterned state (as shown in Fig. 4 in the main text); soon after patterns form in both populations, an oscillatory instability ensues; this linear instability stems from coupling between the spatial modes of the perturbation and those of the steady pattern with matched wavelengths.