Rare event simulation for T-cell activation
The problem of \emph{statistical recognition} is considered, as it arises in immunobiology, namely, the discrimination of foreign antigens against a background of the body's own molecules. The precise mechanism of this foreign-self-distinction, thoug…
Authors: Florian Lipsmeier, Ellen Baake
Noname man usript No. (will b e inserted b y the editor) Florian Lipsmeier · Ellen Baak e Rare ev en t sim ulation for T-ell ativ ation Reeiv ed: date / A epted: date Abstrat The problem of statisti al r e o gnition is onsidered, as it arises in imm unobiology , namely , the disrimination of foreign an tigens against a ba kground of the b o dy's o wn moleules. The preise me hanism of this foreign-self-distintion, though one of the ma jor tasks of the imm une system, on tin ues to b e a fundamen tal puzzle. Reen t progress has b een made b y v an den Berg, Rand, and Burroughs [33 ℄, who mo delled the pr ob abilisti nature of the in teration b et w een the relev an t ell t yp es, namely , T-ells and an tigen-presen ting ells (APCs). Here, the sto hastiit y is due to the random sample of an tigens presen t on the surfae of ev ery APC, and to the random reeptor t yp e that haraterises individual T-ells. It has b een sho wn previously [33 , 37 ℄ that this mo del, though highly idealised, is apable of repro duing imp ortan t asp ets of the reognition phenomenon, and of explaining them on the basis of sto hasti rare ev en ts. These results w ere obtained with the help of a rened large deviation theorem and w ere th us asymptoti in nature. Sim ulations ha v e, so far, b een restrited to the straigh tforw ard simple sampling approa h, whi h do es not allo w for sample sizes large enough to address more detailed questions. Building on the a v ailable large deviation results, w e dev elop an imp ortane sampling te hnique that allo ws for a on v enien t exploration of the relev an t tail ev en ts b y means of sim ulation. With its help, w e in v estigate the me hanism of statistial reognition in some depth. In partiular, w e illustrate ho w a foreign an tigen an stand out against the self ba kground if it is presen t in suien tly man y opies, although no a priori dierene b et w een self and nonself is built in to the mo del. Keyw ords Imm unobiology · statistial reognition · large deviations · rare ev en t sim ulation P A CS 87.16.af · 87.16.dr · 87.18.Tt Mathematis Sub jet Classiation (2000) 92-08 · 92C99 · 60F10 1 In tro dution The notion of statistial reognition b et w een randomly enoun tered moleules is en tral to man y biologial phenomena. This is partiularly eviden t in biologial rep ertoires, whi h on tain enough moleular div ersit y to bind pratially an y randomly enoun tered target moleule. The reep- tor rep ertoire of the imm une system pro vides the b est-kno wn example of a system displa ying probabilit y-based in terations; another one is the olfatory reeptor rep ertoire, whi h reognises m ultitudes of o doran ts. This hane reognition is a w ell-established phenomenon and has b een F. Lipsmeier, E. Baak e F ault y of T e hnology , Bielefeld Univ ersit y , 33501 Bielefeld, German y E-mail: {ipsmei,ebaak e}te hfak.uni-bielefeld.de 2 analysed with the help of v arious statistial and bioph ysial mo dels; ompare [17 , 23 ℄. Here w e will ta kle a mo del of statistial reognition b et w een el l surfa es (in the sense of olletions of n umerous surfae moleules, rather than single ones) of the imm une system. It desrib es a vital prop ert y of our imm une system, whi h omes in to pla y when a virus in v ades the b o dy and starts to m ultiply . F ortunately , ho w ev er, so oner or later it is reognised as a foreign in truder b y ertain white blo o d ells, whi h are part of the imm une system and start a sp ei imm une resp onse that nally eliminates the virus p opulation. This abilit y of the imm une system to disriminate safely b et w een foreign and self moleules is a fundamen tal ingredien t to ev eryda y surviv al of ja w ed v ertebrates; but ho w this w orks exatly is still enigmati. Indeed, the imm une system faes an enormous hallenge b eause it m ust reognise one (or a few) t yp e(s) of (p oten tially dangerous) foreign moleules against an enormous v ariet y of (harmless) self moleules. The partiular diult y lies in the fat that there an b e no a priori dierene b et w een self and nonself (lik e some fundamen tal dierene in moleular struture), sine this w ould op en up the p ossibilit y for moleular mimiry on the part of the pathogen, whi h ould qui kly ev olv e imm uno-in visibilit y b y imitating the self struture. The problem ma y b e phrased as statisti al r e o gnition of one partiular foreign signal against a large, utuating self ba kground. Ho w ev er, imm une biology has b een largely treated deterministially un til, reen tly , an expliit sto hasti mo del w as in tro dued b y v an den Berg, Rand and Burroughs [33 ℄ (heneforth referred to as BRB) and further dev elop ed b y Zin t, Baak e and den Hollander [ 37 ℄. It desrib es (random) enoun ters b et w een the t w o ruial t yp es of white blo o d ells in v olv ed (see Figs. 1 and 2): the an tigen-presen ting ells (APCs), whi h displa y a mixture of self and foreign an tigens at their surfae (a sample of the moleules around in the b o dy), and the T-ells, whi h san the APCs b y means of ertain reeptors and ultimately deide whether or not to reat, i.e., to start an imm une resp onse. APC TCR T cell MHC molecule various peptides internalise digest Fig. 1 A T-ell and an an tigen-presen ting ell (based on Fig. 1 of [35 ℄). An APC absorbs moleules and partiles from its viinit y and breaks them do wn. The emerging fragmen ts, so-alled p eptides (short sequenes of amino aids), serv e as an tigens. They are b ound to so-alled MHC moleules (still within the ell), and the resulting omplexes, ea h omp osed of an MHC moleule and a p eptide, are presen ted on the surfae of the ell (the MHC moleules serv e as arriers or an hors to the ell surfae). Sine most of the moleules in the viinit y of an APC are self moleules, ev ery APC displa ys a large v ariet y of dieren t t yp es of self an tigens and, p ossibly , one (or a small n um b er of ) foreign t yp es. The v arious an tigen t yp es o ur in v arious op y n um b ers. Ea h T- el l is haraterised b y a sp ei t yp e of T-ell reeptor (TCR), whi h is displa y ed in man y identi al opies on the surfae of the partiular T-ell. When a T-ell meets an APC, the on tat b et w een them is established b y a temp orary b ond b et w een the ells, in whi h the TCRs and the MHC-p eptide omplexes in terat with ea h other, whi h results in stim uli to the T-ell b o dy . If the added stim ulation rate is ab o v e a giv en threshold, the T-ell is ativ ated to repro due, and the resulting lones of T-ells will initiate an imm une reation against the in truder. 3 T o b e biologially more preise, w e onsider the enoun ters of so-alled naive T- el ls with pr ofessional APCs in the se ondary lymphoid tissue . A naive T- el l is a ell that has nished its maturation pro ess in the th ym us and has b een released in to the b o dy , where it has not y et b een exp osed to an tigen. It tends to dw ell in se ondary lymphoid tissue lik e lymph no des, where it omes in to on tat with pr ofessional APCs , sp eial white blo o d ells with so-alled MHC moleules at their surfae that serv e as arriers for an tigens. Ea h T- el l is haraterised b y a sp ei t yp e of T-ell reeptor (TCR), whi h is displa y ed in man y identi al opies on the surfae of the partiular T-ell. A large n um b er (estimated at 10 7 in [1℄) of dieren t reeptors, and hene dieren t T-ell t yp es, are presen t in an individual (ev ery t yp e, in turn, is presen t in sev eral opies, whi h form a T-ell lone). Ho w ev er, the n um b er of p oten tial an tigen t yp es is still v astly larger (roughly 10 13 ; see [20 ℄). Th us, sp ei reognition (where one TCR reognises exatly one an tigen) is imp ossible; this is kno wn as Mason's parado x. The task is further ompliated b y the fat that ev ery APC displa ys on the order of thousand(s) of dieren t self an tigen t yp es, in v arious op y n um b ers [ 14 , 20 , 28 ℄, together with, p ossibly , one (or a small n um b er of ) foreign t yp es; the T-ells therefore fae a literal needle in a ha ysta k problem. F or an enoun ter b et w een a pair of T-ell and APC, b oth hosen randomly from the div erse p o ol of T-ells and APCs, the probabilit y to reat m ust b e v ery small (otherwise, imm une reations w ould o ur p ermanen tly); this is a en tral theme in the analysis. It en tails that some questions ma y b e answ ered analytially with the help of large deviation theory; others require sim ulation, but the use of this has b een limited due to the small probabilities in v olv ed, at least with the straigh tforw ard sim ulation metho ds applied so far [ 33 , 37 ℄. The main purp ose of this artile is to devise an eien t imp ortane sampling metho d based on large deviation theory and tailored to the problem at hand, and to use this to in v estigate the me hanism of statistial reognition in more detail. The pap er is organised as follo ws. In Set. 2, w e presen t the most imp ortan t biologial fats and reapitulate the mo del; this will b e a self-on tained, but highly simplied outline, sine the full piture is a v ailable elsewhere [33 , 37 ℄. In Set. 3, w e summarise (mainly from [9℄ and [5℄) some general theory that allo ws to design eien t metho ds of rare ev en t sim ulation on the basis of a large deviation analysis, and tailor these to the problem at hand in Set. 4. Set. 5 presen ts the sim ulation results and analyses them b oth from the omputational and the biologial p oin t of view. Sim ulation sp eeds up b y a fator of nearly 1500 relativ e to the straigh tforw ard approa hes used so far. This enables us to explore regions of parameter spae as y et inaessible, to v alidate previous asymptoti results, and to in v estigate the me hanism of statistial reognition in more depth than previously p ossible. 2 The T-ell mo del In this Setion, w e briey motiv ate and in tro due the mo del of T-ell reognition as rst prop osed b y BRB in 2001 [ 33 ℄ and further dev elop ed b y Zin t, Baak e and den Hollander [ 37 ℄. More preisely , w e only onsider the to y v ersion of this mo del, whi h neglets the mo diation of the T-ell rep ertoire during maturation in the th ym us. This to y v ersion already aptures imp ortan t asp ets of the phenomenon while b eing partiularly transparen t. W e will ome ba k to maturation (already inluded in [33 ℄) in the disussion. When T-ells and APCs meet, the T-ell reeptors bind to the v arious an tigens presen ted b y the APC [6 ℄. F or ev ery single reeptor-an tigen pair, there is an asso iation-disso iation reation, the rate onstan ts for whi h dep end on the mat h of the moleular strutures of reeptor and an tigen. Assuming that asso iation is m u h faster than disso iation and that there is an abundane of reeptors (so that the an tigens are mostly in the b ound state), one an desrib e the reation in terms of the disso iation rates only . Ev ery time a reeptor un binds from an an tigen, it sends a signal to the T-ell, pro vided the asso iation has lasted for at least one time unit (i.e., w e resale time so that the unit of time is this minimal asso iation time required). The duration of a binding of a giv en reeptor-an tigen pair follo ws the Exp (1 /τ ) distribution, i.e. the exp onen tial distribution with mean τ , where τ is the in v erse disso iation rate of the pair in question. The rate of stim uli indued b y the in teration of 4 APC 2 T−cell 2 T−cell 2 APC 1 T−cell 1 APC 3 T−cell 3 Fig. 2 Cariature of T-ells and APCs (from [37 ℄). Note that ev ery T-ell has man y opies of one partiular reeptor t yp e, but dieren t T-ells ha v e dieren t reeptor t yp es. In on trast, ev ery APC arries a mixture of an tigen t yp es, whi h ma y app ear in v arious op y n um b ers. our an tigen with the reeptors in its viinit y is then giv en b y w ( τ ) = 1 τ exp( − 1 τ ) , (1) i.e., the disso iation rate times the probabilit y that the asso iation has lasted long enough. (If the simplifying assumption of unlimited reeptor abundane is disp ensed with, Eq. ( 1) m ust b e mo died, see [34℄.) As sho wn in Fig. 3 , the funtion w rst inreases and then dereases with τ with a maxim um at τ = 1 , whi h reets the fat that, for τ < 1 , the bindings tend not to last long enough, whereas for τ > 1 , they tend to last so long that only few stim uli are exp eted p er time unit. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 2 4 6 8 10 w ( τ ) τ 0 2 4 6 8 10 12 14 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 densit y of W densit y of W ϑ Fig. 3 Left: the funtion w . Righ t: the densities of W = w ( T ) and W ϑ with tilting parameter ϑ = 46 (f. Set. 3.2 ). The densities ha v e p oles at w ( 0) = 0 and w ( 1) = 0 . 367 9 (due to the v anishing deriv ativ e of w at τ = 0 and τ = 1 ), but the righ t p oles are in visible b eause they supp ort v ery little probabilit y mass. In fat, for ε = 0 . 01 , one has P (0 ≤ W ≤ ε ) = 0 . 98 and P ( w (1) − ε ≤ W ≤ w (1)) = 2 . 17 · 10 − 9 , whereas P (0 ≤ W ϑ ≤ ε ) = 0 . 137138 and P ( w (1) − ε ≤ W ϑ ≤ w (1)) = 0 . 0050 . The T-ell sums up the signals indued b y the dieren t an tigens on the APC, and if the total stim ulation rate rea hes a ertain threshold v alue, the ell initiates an imm une resp onse. This mo del relies on sev eral h yp otheses, whi h are kno wn as kineti pro ofreading [21 , 22, 18 , 13 ℄, serial triggering [31 , 30 , 27, 4, 29 , 10 ℄, oun ting of stim ulated TCRs [ 36 , 25 ℄, and the optimal dw ell-time h yp othesis [15 , 12 ℄. 5 Due to the h uge amoun t of dieren t reeptor and an tigen t yp es, it is imp ossible (and unne- essary) to presrib e the binding durations for all pairs of reeptor and an tigen t yp es individually . Therefore, BRB hose a probabilisti approa h to desrib e the meeting of APCs and T-ells. A randomly hosen T-ell (that is, a randomly hosen t yp e of reeptor) enoun ters a randomly ho- sen APC (that is, a random mixture of an tigens). The mean binding time that go v erns the binding of this random reeptor to the j th t yp e of an tigen is tak en to b e a random v ariable denoted b y T j . The T j are indep enden t and iden tially distributed (i.i.d.) and are assumed to follo w the Exp(1 / ¯ τ ) distribution, i.e., the exp onen tial distribution with mean ¯ τ , where ¯ τ is a free parameter. Note that there are t w o exp onen tial distributions (and t w o lev els of a v eraging) in v olv ed here. First, the duration of an individual binding b et w een a t yp e- j an tigen and a random reeptor is Exp(1 / T j ) distributed (see the disussion of Eq. (1)). Seond, T j , the me an dur ation of su h a binding (where the reeptor is hosen one and the times are a v eraged o v er rep eated bindings with a j an tigen) is itself an exp onen tial random v ariable, with realisation τ j . Finally , its mean, E ( T j ) = ¯ τ , is the mean binding time of a j -an tigen (and, due to the i.i.d. assumption, of an y an tigen) when a v eraged o v er all enoun ters with the v arious reeptor t yp es. The exp onen tial distribution of the individual binding time is an immediate onsequene of the (rst-order) un binding kinetis. In on trast, the orresp onding assumption for the T j is made for simpliit y; the approa h is ompatible with v ar- ious other distributions as w ell, see [33 ℄ and [ 37 ℄. The i.i.d. assumption, ho w ev er, is ruial, sine it implies, in partiular, that there is no dierene b et w een self and foreign an tigens here; i.e., no a priori distintion is built in to the mo del. The total stim ulation a T-ell reeiv es is the sum o v er all stim ulus rates W j = w ( T j ) that emerge from an tigens of the j 'th t yp e. It is further assumed that there is at most one t yp e of foreign an tigen in z ( f ) opies on an APC, whose signal m ust b e disriminated against the signals of a h uge amoun t of self an tigens. (There ould, in priniple, b e m ultiple foreign p eptide t yp es, but there are go o d reasons to assume that there are me hanisms to ensure that a giv en T-ell sees at most one foreign p eptide t yp e, see [34 ℄). The self an tigens are here divided in to t w o distint lasses, c and v , that are presen t in dieren t op y n um b ers z ( c ) and z ( v ) . An APC displa ys m ( c ) and m ( v ) dieren t t yp es of lass c and v . The indies c and v stand for onstitutiv e and for v ariable, resp etiv ely; but for the purp ose of this artile, only the abundanies are relev an t, in partiular, z ( c ) > z ( v ) and m ( c ) < m ( v ) . Ov er the whole APC the total n um b er of an tigens is then m ( c ) z ( c ) + m ( v ) z ( v ) =: M if no foreign an tigen is presen t. If z ( f ) foreign moleules are also presen t, the self moleules are assumed to b e prop ortionally displaed (via the fator q := ( M − z ( f ) ) / M ), so that the total n um b er of an tigens remains un hanged at z ( f ) + m ( c ) q z ( c ) + m ( v ) q z ( v ) = M . (2) The total stim ulation rate in a random enoun ter of T-ell and APC an then b e desrib ed as a funtion of z ( f ) : G ( z ( f ) ) := m ( c ) X j =1 q z ( c ) W j + m ( c ) + m ( v ) X j = m ( c ) +1 q z ( v ) W j + z ( f ) W m ( c ) + m ( v ) +1 , (3) i.e., a w eigh ted sum of i.i.d. random v ariables. Alternativ ely , w e onsider the extension of the mo del prop osed b y Zin t et al. [37 ℄, whi h, instead of the deterministi op y n um b ers z ( c ) , z ( v ) , uses random v ariables Z ( c ) j , Z ( v ) j distributed aording to binomial distributions with E ( Z ( c ) j ) = z ( c ) , E ( Z ( v ) j ) = z ( v ) , where E denotes exp etation (so the exp eted n um b er of an tigens p er APC is still M ). The mo del then reads G ( z ( f ) ) := m ( c ) X j =1 q Z ( c ) j W j + m ( c ) + m ( v ) X j = m ( c ) +1 q Z ( v ) j W j + z ( f ) W m ( c ) + m ( v ) +1 . (4) In line with [33 , 37 ℄, w e n umerially sp eify the mo del parameters as follo ws: ¯ τ = 0 . 04 ; m ( c ) = 50 , m ( v ) = 1500 , z ( c ) = 500 , z ( v ) = 50 (and hene M = 10 5 ). The distributions in the extended 6 mo del are the binomials Bin( ζ ( c ) , p ) and Bin( ζ ( v ) , p ) for Z ( c ) j and Z ( v ) j resp etiv ely , where ζ ( c ) = 1000 , ζ ( v ) = 100 , a nd p = 0 . 5 . The relev an t quan tit y for us is no w the probabilit y P G ( z ( f ) ) ≥ g act (5) that the stim ulation rate rea hes or surpasses a threshold g act . T o a hiev e a go o d foreign-self disrimination, there m ust b e a large dierene in probabilit y b et w een the stim ulation rate in the ase with self an tigens only ( z ( f ) = 0 ), and the stim ulation rate with the foreign an tigen presen t, i.e., 1 ≫ P G ( z ( f ) ) ≥ g act ≫ P G (0) ≥ g act ≥ 0 (6) for realisti v alues of z ( f ) . Note that b oth ev en ts m ust b e rare ev en ts otherwise, the imm une system w ould re all the time. Th us g act m ust b e m u h larger than E ( G ( z ( f ) )) (whi h, due to (2 ) and the iden tial distribution of the W j , is indep enden t of z ( f ) ). Ev aluating these small probabilities is a hallenge. So far, t w o routes ha v e b een used: analyti (asymptoti) theory based on large deviations (LD) and straigh tforw ard sim ulation (so-alled simple sampling). Both ha v e their shortomings: the LD approa h is only exat in the limit of innitely man y an tigen t yp es (and the a v ailable error estimates are usually to o rude to b e useful); the sim ulation strategy , on the other hand, is so time-onsuming that it b eomes simply imp ossible to obtain sample sizes large enough for a detailed analysis, in partiular for large v alues of g act . Therefore, an imp ortane sampling approa h is required. Let us no w reapitulate some underlying theory . 3 Rare ev en t sim ulation: general theory The general problem w e no w onsider is to estimate the probabilit y P ( A ) of a (rare) ev en t A under a probabilit y measure P . The straigh tforw ard approa h, kno wn as simple sampling, uses the estimate ( [ P ( A )) N := 1 N N X i =1 1 { S ( i ) ∈ A } = 1 N card { 1 ≤ i ≤ N | S ( i ) ∈ A } , (7) where the { S ( i ) } 1 ≤ i ≤ N are indep enden t and iden tially distributed (i.i.d.) random v ariables with distribution P , 1 { . } denotes the indiator funtion, and N is the sample size; w e will throughout use b v for an estimate of a quan tit y v . ( [ P ( A )) N is ob viously an un biased and onsisten t estimate, but, for small P ( A ) , the on v ergene to P ( A ) is slo w, and large samples are required to get reliable estimates. V arious sim ulation metho ds are a v ailable that deal with this problem and yield a b etter rate of on v ergene (see the monograph b y Bu klew [5℄ for an o v erview). Most of them a hiev e this impro v emen t b y reduing the v ariane of the estimator. W e will onen trate here on the most wide-spread lass of metho ds, namely imp ortane sampling. As is w ell kno wn, one in tro dues a new sampling distribution Q here under whi h A is more lik ely to happ en, pro dues samples from this distribution and returns to the original distribution b y rew eigh ting. In general, nding a go o d imp ortane sampling distribution that redues the v ariane as m u h as p ossible is an art, and m u h of the literature rev olv es around this. Some general purp ose and man y ad ho strategies exist, but usually , imp ortane sampling distributions are b est tailored b y exploiting the struture of the sp ei problem at hand. Ho w ev er, if the problem an b e em b edded in to a sequene of problems for whi h a so-alled large deviation priniple is v alid, a unied theory is a v ailable that iden ties the most eien t sim ulation distribution. This te hnique of large deviation sim ulation w as in tro dued b y Sado wski and Bu klew [ 26 ℄, laid do wn in the monograph b y Bu klew [5 ℄, and further dev elop ed b y Diek er and Mandjes [ 9℄. It rests on the w ell-established theory of large deviations, as summarised, for example, in the b o oks b y Dem b o and Zeitouni [ 7℄ or den Hollander [8℄. Let us reapitulate the basi ba kground. 7 3.1 Large deviation probabilities Consider a sequene { S n } of random v ariables on the probabilit y spae ( R d , B , P ) , where B is the Borel σ -algebra of R d . Let { P n } b e the family of probabilit y measures indued b y { S n } , i.e., P n ( B ) = P ( S n ∈ B ) for B ∈ B . W e assume throughout that { S n } satises a large deviation priniple (LDP) aording to the follo wing denition [7 , 9℄: Denition 1 (Large deviation priniple) A family of probabilit y measures { P n } on ( R d , B ) satises the large deviation priniple (LDP) with rate funtion I if I : R d → [0 , ∞ ] is lo w er semion tin uous and, for all B ∈ B , − inf x ∈ B ◦ I ( x ) ≤ lim inf n →∞ 1 n log P n ( B ) ≤ lim sup n →∞ 1 n log P n ( B ) ≤ − inf x ∈ B I ( x ) , (8) where B ◦ := int( B ) and B := clos( B ) denote the in terior and the losure of B , resp etiv ely . I is said to b e a go o d rate funtion if it has ompat lev el sets in that I − 1 ([0 , c ]) = { x ∈ R d : I ( x ) ≤ c } is ompat for all c ∈ R d . ⊓ ⊔ A set B is alled an I - ontinuity set if inf x ∈ B ◦ I ( x ) = inf x ∈ B I ( x ) = inf x ∈ B I ( x ) . (9) If B is su h a set, the LDP means that P n ( B ) dea ys exp onen tially for large n , with dea y o eien t inf x ∈ B I ( x ) . A p oin t b is alled a minimum r ate p oint of B if inf x ∈ B I ( x ) = I ( b ) . Large deviation priniples are w ell kno wn for man y families of random v ariables, lik e empirial means of i.i.d. random v ariables or empirial measures of Mark o v hains. F or the appliation w e ha v e in mind, whi h in v olv es sums of indep enden t, but not iden tially distributed random v ariables, w e need the fairly general setting of the Gärtner-Ellis theorem, whi h w e reapitulate here (f. [7, Thm. 2.3.6℄ and [8 , Ch. V℄). Let ϕ n ( ϑ ) := E P n ( e h ϑ,S n i ) , ϑ ∈ R d , b e the momen t-generating funtion of S n , where h ., . i denotes the salar pro dut and E µ ( . ) denotes the exp etation of a random v ariable with resp et to the probabilit y measure µ . Theorem 1 (Gärtner-Ellis) Assume that (G1) lim n →∞ 1 n log ϕ n ( nϑ ) =: Λ ( ϑ ) ∈ [ −∞ , ∞ ] exists , (G2) 0 ∈ int ( D Λ ) , wher e D Λ := { ϑ ∈ R d : Λ ( ϑ ) < ∞} is the eetiv e domain of Λ , (G3) Λ is lower semi- ontinuous on R d , (G4) Λ is dier entiable on int ( D Λ ) , (G5) Either D Λ = R d or Λ is ste ep at its b oundary ∂ D Λ , i.e., lim int( D Λ ) ∋ ϑ → ∂ D Λ |∇ Λ ( ϑ ) | = ∞ . Then, { P n } satises the LDP on R d with go o d r ate funtion I , wher e I is the L e gendr e tr ansform of Λ , i.e., I ( x ) = sup ϑ ∈ R d [ h x, ϑ i − Λ ( ϑ )] , x ∈ R d . (10) ⊓ ⊔ The funtion Λ in (G1) is on v ex. If there is a solution ϑ ∗ of ∇ Λ ( ϑ ) = x, (11) one has I ( x ) = h ϑ ∗ , x i − Λ ( ϑ ∗ ) . (12) If Λ is stritly on v ex in all diretions, ϑ ∗ is unique. See Fig. 4 for a one-dimensional example (the T-ell appliation, in fat). 8 3.2 Sim ulating rare ev en t probabilities Let no w A ∈ B b e a r ar e event in the sense that 0 < inf x ∈ A I ( x ) < ∞ . Here, the rst inequalit y implies that A b eomes exp onen tially unlik ely as n → ∞ , whereas the seond inequalit y serv es to exlude nongeneri ases (in partiular ases where the ev en t is imp ossible). An imp ortan t notion for the rare ev en t sim ulation of P n ( A ) is that of a dominating p oint [ 5, p. 83℄: A p oin t a is a dominating p oint of the set A if it is the unique p oin t su h that a) a ∈ ∂ A , b) ∃ a unique solution ϑ ∗ of ∇ Λ ( ϑ ) = a , and ) A ⊂ { x ∈ R d : h ϑ ∗ , x − a i ≥ 0 } . A dominating p oin t, if it exists, is alw a ys a unique minim um rate p oin t (see [ 5, p. 83℄). Con- v exit y of A implies existene of a dominating p oin t (f. [9 ℄). F ollo wing [9℄ w e no w turn to the problem of sim ulating P n ( A ) = E P n ( 1 { S n ∈ A } ) . The naiv e simple-sampling estimate obtained from N i.i.d. opies S ( i ) n ( 1 ≤ i ≤ N ), dra wn from P n , is, as in (7), giv en b y \ P n ( A ) N := 1 N N X i =1 1 { S ( i ) n ∈ A } . (13) It is un biased and on v erges (almost surely) to P n ( A ) in the limit N → ∞ , but it is ineien t sine it requires that N inrease exp onen tially with n to yield a meaningful estimate. Instead of { S n } , one therefore onsiders an alternativ e family of random v ariables, { T n } with distribution family { Q n } , again on ( R d , B ) , under whi h A o urs more frequen tly . Assuming that P n and Q n are absolutely on tin uous with resp et to ea h other, one an use the iden tit y P n ( A ) = E P n ( 1 { S n ∈ A } ) = E Q n 1 { T n ∈ A } d P n d Q n ( T n ) , (14) where d P n / d Q n is the Radon-Nik o dym deriv ativ e of P n with resp et to Q n . The resulting imp or- tane sampling estimate then relies on i.i.d. samples T ( i ) n from { Q n } and reads \ P Q n ( A ) N := 1 N N X i =1 1 { T ( i ) n ∈ A } d P n d Q n ( T ( i ) n ) , (15) where (d P n / d Q n )( . ) ats as a rew eigh ting fator from the sampling distribution to the original one. It is reasonable to assume that (d P n / d Q n ) is on tin uous to a v oid the usual problems with L 1 -funtions; this is no restrition for our in tended appliation. An adequate optimalit y onept in this on text is that of asymptoti eieny . A ording to [9℄, it is based on the r elative err or η N ( Q n , A ) dened via its square η 2 N ( Q n , A ) := V Q n \ P Q n ( A ) N P n ( A ) 2 (16) (where V µ ( . ) denotes the v ariane of a random v ariable with resp et to the probabilit y measure µ ). The relativ e error is prop ortional to the width of the ondene in terv al relativ e to the (exp eted) estimate itself. Asymptoti eieny is then dened as follo ws. Denition 2 (Asymptoti eieny) An imp ortane sampling family { Q n } is alled asymp- toti al ly eient for the rare ev en t A if lim n →∞ 1 n log N ∗ Q n = 0 , (17) where N ∗ Q n := inf { N ∈ N : η N ( Q n , A ) ≤ η max } for some giv en maximal relativ e error η max , 0 < η max < ∞ . 9 In w ords, asymptoti eieny means that the n um b er of samples required to k eep the relativ e error b elo w a presrib ed b ound η max inreases only sub exp onen tially (rather than exp onen tially as with simple sampling). The onrete hoie of η max is atually irrelev an t, see Lemma 1 in [9℄. An ob vious idea from large deviation theory w ould b e to use, as sampling distributions, the family of measures { P ϑ n } that are exp onen tially tilted with parameter ϑ , that is, d P ϑ n d P n ( x ) = e n h ϑ,x i ϕ n ( nϑ ) , x ∈ R d ; (18) P ϑ n then tak es the role of Q n . The task remains to nd a suitable ϑ , i.e., a tilting parameter that mak es { P ϑ n } asymptotially eien t. Neessary and suien t onditions for this are giv en in [9, Assumption 1 and Corollary 1℄ and are summarised b elo w, in a form adapted to the presen t on text. Theorem 2 (Diek er-Mandjes 2005) Assume that, for some given ϑ ∗ , (V1) { P n } satises an LDP with go o d r ate funtion I , (V2) lim sup n →∞ 1 n log ϕ n ( γ nϑ ∗ ) < ∞ for some γ > 1 , and, likewise, with ϑ ∗ r epla e d by − ϑ ∗ , (V3) The r ar e event A is b oth an I - ontinuity set and an ( I + h ϑ ∗ , . i ) - ontinuity set. Then, the tilte d me asur e { P ϑ ∗ n } is asymptoti al ly eient for simulating A if and only if inf x ∈ R d [ I ( x ) − h ϑ ∗ , x i ] + inf x ∈ A [ I ( x ) + h ϑ ∗ , x i ] = 2 inf x ∈ A ◦ I ( x ) . (19) W e use assumption (V2) here to replae the w eak er but less easy to v erify ondition (2) in Assumption 1 of [ 9 ℄, in line with the paragraph b elo w (2) in [ 9℄, or [7, Thm. 4.3.1℄. Note also that (V2) holds automatially if ϕ n ( nϑ ) exists for all ϑ but this is not mandatory here, sine only a giv en ϑ ∗ is onsidered. The pro of of Theorem 2 is giv en in [9℄ and need not b e reapitulated here; but w e w ould lik e to ommen t briey on what happ ens in the en tral ondition (19 ). Replaing Q n b y P ϑ ∗ n in (16) and (15), w e an rewrite η 2 N as η 2 N ( P ϑ ∗ n , A ) = V P ϑ ∗ n ( \ P P ϑ ∗ n ( A )) N P n ( A ) 2 = 1 N V P ϑ ∗ n ( \ P P ϑ ∗ n ( A )) 1 P n ( A ) 2 = 1 N 1 P n ( A ) 2 h Z A d P n d P ϑ ∗ n 2 d P ϑ ∗ n − P n ( A ) 2 i . (20) Ob viously (b y (V1) and (V3)), 2 inf x ∈ A ◦ I ( x ) (i.e., the righ t-hand side of (19 )) is the exp onen tial dea y rate of ( P n ( A )) 2 . Insp etion of the pro of of Theorem 2 rev eals that the left-hand side of (19 ) is the exp onen tial dea y rate of R A d P n d P ϑ ∗ n 2 d P ϑ ∗ n . It is lear from (20 ) that, for asymptoti eieny to hold, R A d P n d P ϑ ∗ n 2 d P ϑ ∗ n m ust tend to zero at least as fast as ( P n ( A )) 2 . But it annot derease faster, sine V P ϑ ∗ n ( \ P P ϑ ∗ n ( A )) 1 is nonnegativ e, so that R A d P n d Q n 2 d Q n ≥ ( P n ( A )) 2 for arbitrary Q n . Hene, the exp onen tial dea y rates m ust b e exatly equal, as stated b y (19 ). (A losely related argumen t is giv en in [5, Ch. 5.2℄.) Theorem 2 is widely appliable. It holds in man y standard situations, in partiular in man y of those that arise in appliations. Prop osition 1 L et { P n } b e a family of pr ob ability me asur es that satisfy the onditions of the Gärtner-El lis the or em, with (go o d) r ate funtion I . L et A b e a r ar e event with dominating p oint a , let ϑ ∗ b e the unique solution of ∇ Λ ( ϑ ) = a , and assume (V2) and (V3) . Then { P ϑ ∗ n } is the unique tilte d family that is asymptoti al ly eient for simulating P n ( A ) . 10 Pr o of The pro of is a simple appliation of Thm. 2. (V1) follo ws from the Gärtner-Ellis theorem; w e only need to v erify ondition (19 ). F or the rst inm um in (19 ), one obtains inf x ∈ R d [ I ( x ) − h ϑ ∗ , x i ] = − Λ ( ϑ ∗ ) = I ( a ) − h ϑ ∗ , a i . (21) Here, the rst step follo ws from the on v ex dualit y lemma (ompare [7 , Lemma 4.5.8℄), whi h is appliable sine Λ is lo w er semion tin uous b y (G3), and on v ex and > −∞ ev erywhere (this follo ws from (G1) and (G2) b y [ 8, Lemma V.4℄). The seond step is due to part b) of the dominating p oin t prop ert y of a , together with Eq. ( 12 ). As to the seond inm um in (19 ), a minimises b oth I and h ϑ ∗ , . i on A (b y the dominating p oin t prop ert y). T ogether with (V3), this giv es inf x ∈ A [ I ( x ) + h ϑ ∗ , x i ] = I ( a ) + h ϑ ∗ , a i . (22) Eqs. (21 ) and (22) together giv e (19 ) b eause inf x ∈ A ◦ I ( x ) = inf x ∈ ∂ A I ( x ) = I ( a ) . ⊓ ⊔ R emark 1 Note that an eieny result losely related to Prop osition 1 has previously b een giv en b y Bu klew [5 , Thm. 5.2.1℄, but this is based on the v ariane rather than the relativ e error; and it is only a suien t ondition. Note also that our assumption of a dominating p oin t greatly simplies the situation. Theorem 2 also allo ws to op e with situations without a dominating p oin t but this is not needed b elo w. Let us no w apply this theory to the T-ell mo del. 4 Rare ev en t sim ulation: the T-ell mo del Reall that sim ulating the T-ell mo del means sampling the random v ariables G ( z ( f ) ) of (3) and estimating the orresp onding tail probabilities P ( G ( z ( f ) ) ≥ g act ) . Insp etion of Eq. ( 3 ) rev eals t w o diulties: 1. G ( z ( f ) ) is a w eigh ted sum of i.i.d. random v ariables, to whi h the standard results for sums of i.i.d. random v ariables (in partiular, Cramér's theorem) are not appliable. W e therefore need an extension to w eigh ted sums or, b etter, to general sums of indep enden t, but not iden tially distributed random v ariables, whi h inlude w eigh ted sums as a simple sp eial ase. This is straigh tforw ard and will b e the sub jet of Set. 4.1 . In partiular, it will b e seen that, lik e in the i.i.d. ase, ev ery term in the sum m ust b e tilted with the same parameter, but no w this global tilting fator is a funtion of all the individual distributions in v olv ed. 2. Sim ulating the random v ariables W j = w ( T j ) is straigh tforw ard via simple sampling: dra w Exp(1 / ¯ τ ) distributed random n um b ers τ j (as realisations of T j ) and apply the transformation (1 ) . Ho w ev er, sim ulating the orresp onding tilted v ariables is a diult task, for t w o reasons. First of all, there is no indiation of ho w to sample from the tilted distribution via transfor- mation of one of the elemen tary distributions (lik e Uni [0 , 1] (the uniform distribution on the unit in terv al), or Exp( λ )) for whi h eien t random n um b er generation is p ossible. Although su h a transformation migh t exist in priniple, there is no systemati w a y of nding it. One reason for this is that tilting ats at the lev el of the densities, but ev en the original (un tilted) densit y of W = w ( T ) is not a v ailable expliitly . (With W and T (without indies) w e mean an y represen tativ e of the family .) This is b eause its alulation requires the in v erse funtions and deriv ativ es of the t w o bran hes (inreasing and dereasing) of the funtion w , but these are una v ailable analytially . In the absene of a transformation metho d, one migh t onsider to determine the tilted densit y n umerially , in tegrate it (again n umerially) and disretise and tabulate the resulting distri- bution funtion. Ho w ev er, this is, again, forbidding for our partiular funtion w : due to the v anishing deriv ativ es at T = 0 and T = 1 , the transformation form ula for densities yields singularities in the densit y of W at these v alues, with a sizeable fration of the probabilit y 11 mass onen trated v ery lose to 0 (see Fig. 3 ). This renders n umerial alulations unreliable. T o irum v en t these problems, w e will, in Set. 4.2 , presen t a sampling metho d for the tilted random v ariable W ϑ that is based on tilting T rather than W itself. 4.1 Large deviations for indep enden t but not iden tially distributed random v ariables W e onsider K indep enden t families of i.i.d. R d -v alued random v ariables, { Y (1) ℓ } , . . . , { Y ( K ) ℓ } (i.e., the distribution within an y giv en family { Y ( k ) ℓ } , 1 ≤ k ≤ K , is xed, but the distributions ma y v ary aross families). Assume that Λ ( k ) ( ϑ ) := log E ( e h ϑ,Y ( k ) 1 i ) , the log momen t-generating funtion of Y ( k ) 1 , is nite for all ϑ ∈ R d and 1 ≤ k ≤ K (here, E ( . ) refers to the probabilit y measure indued b y the random v ariable in v olv ed). Let n (1) , . . . , n ( K ) b e p ositiv e in tegers, n := P K k =1 n ( k ) , V n := n (1) X ℓ =1 Y (1) ℓ + . . . + n ( K ) X ℓ =1 Y ( K ) ℓ , (23) and P n b e the probabilit y measure indued b y S n = V n /n . In the limit n → ∞ , sub jet to n ( k ) /n → γ ( k ) for all 1 ≤ k ≤ K , the limiting log-momen t generating funtion of { S n } b eomes Λ ( ϑ ) = lim n →∞ 1 n log E ( e h ϑ,V n i ) = lim n →∞ K X k =1 n ( k ) n Λ ( k ) ( ϑ ) = K X k =1 γ ( k ) Λ ( k ) ( ϑ ) , (24) where the seond step is due to indep endene. Sine, b y assumption, Λ ( k ) ( ϑ ) < ∞ for all ϑ ∈ R d and 1 ≤ k ≤ K , the Λ ( k ) are dieren tiable on all of R d (see [7 , Lemma 2.2.31℄); in fat, they are ev en C ∞ ( R d ) [7 , Ex.erise 2.2.24℄. Th us, Λ is C ∞ ( R d ) as w ell. By (24), w e ha v e (G1). Again due to Λ ( k ) ( ϑ ) < ∞ , (G2) and (G5) are automatially satised. F urthermore, the dieren tiabilit y of Λ en tails (G3) and (G4). W e ha v e therefore sho wn Lemma 1 Under the assumptions of this p ar agr aph, { P n } satises the Gärtner-El lis the or em, with r ate funtion I given by Eq. (10) . ⊓ ⊔ Su h { P n } are therefore andidates for eien t sim ulation aording to Prop. 1 . The tilting fator ϑ ∗ ma y not b e aessible analytially , but an b e ev aluated n umerially from ( 11 ) . Due to inde- p endene, tilting of S n with nϑ ∗ (that is, tilting of V n with ϑ ∗ ) is equiv alen t to tilting ea h Y ( k ) ℓ with ϑ ∗ . 4.2 Tilting of transformed random v ariables Unlik e the W j , the Exp(1 / ¯ τ ) -distributed random v ariables T j are tilted easily (tilting with ϑ simply giv es Exp( − ϑ + 1 / ¯ τ ) ). One is therefore tempted to tilt the T j rather than the W j , or, in other w ords, to in ter hange the order of tilting and transformation. The follo wing Theorem states the k ey idea. Theorem 3 L et X b e an R d -value d r andom variable with pr ob ability me asur e µ , and let Y := h ◦ X (or Y = h ( X ) by slight abuse of notation), wher e h : R d → R d is µ -me asur able. Then Y has pr ob ability me asur e ν = µ ◦ h − 1 , wher e h − 1 ( y ) denotes the pr eimage of y . Assume now that E µ ( e h ϑ,h ( X ) i ) exists, let ˜ X ϑ b e an R d -value d r andom variable with pr ob ability me asur e ˜ µ ϑ r elate d to µ via d ˜ µ ϑ d µ ( x ) = e h ϑ,h ( x ) i E µ ( e h ϑ,h ( X ) i ) (25) 12 (so that ˜ µ ϑ ≪ µ ), and let ˜ Y ϑ = h ( ˜ X ϑ ) . Then, the me asur es ˜ ν ϑ (of ˜ Y ϑ ) and ν ϑ (for the tilte d version of ν , b elonging to Y ϑ ) ar e e qual, wher e ν ϑ ≪ ν with R adon-Niko dym density d ν ϑ d ν ( y ) = e h ϑ,y i E ν ( e h ϑ,Y i ) . (26) Pr o of Note rst that e h ϑ,y i is learly µ -measurable, and E ν ( e h ϑ,Y i ) = Z R d e h ϑ,y i d ν ( y ) = Z R d e h ϑ,h ( x ) i d µ ( x ) = E µ ( e h ϑ,h ( X ) i ) , (27) whi h exists b y assumption, so ν ϑ is w ell-dened. W e no w ha v e to sho w that ˜ ν ϑ ( B ) = ν ϑ ( B ) for arbitrary Borel sets B . Observing that ˜ ν ϑ = ˜ µ ϑ ◦ h − 1 and emplo ying the form ulas for transformation of measures [3 , (13.7)℄ and hange of v ariable [3, Thm. 16.13℄, together with ( 25 ) , one indeed obtains ˜ ν ϑ ( B ) = ˜ µ ϑ h − 1 ( B ) = Z h − 1 ( B ) d ˜ µ ϑ d µ ( x )d µ ( x ) = 1 E µ ( e h ϑ,h ( X ) i ) Z h − 1 ( B ) e h ϑ,h ( x ) i d µ ( x ) = 1 E ν ( e h ϑ,Y i ) Z B e h ϑ,y i d ν ( y ) = Z B d ν ϑ d ν ( y )d ν ( y ) = ν ϑ ( B ) , (28) whi h pro v es the laim. ⊓ ⊔ In w ords, Theorem 3 is nothing but the simple observ ation that, to obtain the tilted v ersion of Y = h ( X ) , one an rew eigh t the measure µ of X with the fators e h ϑ,h ( x ) i , rather than rew eigh ting the measure ν of Y with e h ϑ,y i . It should b e lear, ho w ev er, that the measure ˜ µ ϑ diers from the usual tilted v ersion of µ , whi h w ould in v olv e tilting fators e h ϑ,x i rather than e h ϑ,h ( x ) i ; for this reason, w e use the notation ˜ µ ϑ rather than µ ϑ . Su h kind of tilting is ommon in large deviation theory (see, e.g., [7 , Chap. 2.1.2℄). Nev ertheless, the simple observ ation ab o v e is the k ey to sim ulation if µ (and ˜ µ ϑ ) are readily aessible at least n umerially , but ν (and ν ϑ ) are not. This is preisely our situation, with ˜ T ϑ , αW ϑ and αw ( α ∈ { q z ( c ) , q z ( v ) , z ( f ) } ), resp etiv ely , taking the roles of ˜ X ϑ , Y ϑ and h (w e will use f , ˜ f ϑ , g and g ϑ for the orresp onding densities of T , ˜ T ϑ , αW , and ( αW ) ϑ ). Still, rew eigh ting of the exp onen tial densit y of T with e ϑαw ( τ ) do es not yield an expliit losed-form densit y , and no diret sim ulation metho d is a v ailable for the orresp onding random v ariables. Ho w ev er, the rew eigh ted densities are easily aessible n umerially , in on trast to those of W and its tilted v arian t, W ϑ . The problem ma y th us b e solv ed b y alulating and in tegrating ˜ f ϑ n umerially and disretising and tabulating the resulting distribution funtion ˜ F ϑ . Samples of ˜ T ϑ ma y then b e dra wn aording to this table (i.e., b y formally lo oking up the solution of ˜ F ( ˜ T ϑ ) = U for U ∼ Uni [0 , 1] ), and αW ϑ = αw ( ˜ T ϑ ) is then readily ev aluated. The only diult y left is the time required for sear hing the table. But this is a pratial matter and will b e dealt with in the next paragraph. 4.3 The algorithm T aking together our theoretial results, w e an no w detail the sp ei imp ortane sampling algo- rithm for the sim ulation of the T-ell mo del of Set. 2 . If not stated otherwise, w e will refer to the basi mo del (3). Reall that it desrib es the stim ulation rate G ( z ( f ) ) and w e wish to ev aluate the probabilit y P ( G ( z ( f ) ) ≥ g act ) . T o apply LD sampling, let us em b ed the mo del in to a sequene of mo dels with inreasing total n um b er n = n ( c ) + n ( v ) + n ( f ) of an tigen t yp es, where n ( c ) , n ( v ) , and n ( f ) are the n um b ers of onstitutiv e, v ariable and foreign an tigen t yp es. (This is an aritiial sequene of mo dels required 13 to form ulate the limiting pro ess in v olv ed in the theory; in on trast to the original mo del, there an no w b e m ultiple foreign an tigen t yp es.) Let G n ( z ( f ) ) = n ( c ) X j =1 q n z ( c ) W j + n ( c ) + n ( v ) X j = n ( c ) +1 q n z ( v ) W j + n ( c ) + n ( v ) + n ( f ) X j = n ( c ) + n ( v ) +1 z ( f ) W n ( c ) + n ( v ) + j , (29) where q n = n ( c ) z ( c ) + n ( v ) z ( v ) − n ( f ) z ( f ) n ( c ) z ( c ) + n ( v ) z ( v ) (30) (where z ( c ) , z ( v ) , and z ( f ) are indep enden t of n ). Clearly , G n ( z ( f ) ) oinides with G ( z ( f ) ) of (3) if n ( c ) = m ( c ) , n ( v ) = m ( v ) , and n ( f ) = m ( f ) , where m ( f ) = 0 or m ( f ) = 1 dep ending on whether z ( f ) = 0 or z ( f ) > 0 ; then, n = m = m ( c ) + m ( v ) + m ( f ) . W e ha v e to onsider P G n ( z ( f ) ) /n > g act /m (this reets the fat that g act m ust sale with system size). The sequenes { G n ( z ( f ) ) } and { G n ( z ( f ) ) } /n tak e the roles of { V n } and { S n } , resp etiv ely , in Ses. 3.1 and 4.1 , with P n the la w of G n ( z ( f ) ) /n ; and w e onsider A = [ g act /m, ∞ ) with E ( G m ( z ( f ) ) /m ) < g act /m < M w (1) /m (the latter is the maxim um v alue of G m ( z ( f ) ) /m sine w ( τ ) has its maxim um at τ = 1 ). The limit n → ∞ is then tak en so that lim n →∞ n ( c ) /n = m ( c ) /m , lim n →∞ n ( v ) /n = m ( v ) /m , as w ell as lim n →∞ n ( f ) /n = m ( f ) /m , that is, the relativ e amoun ts of onstitutiv e, v ariable, and foreign an tigens approa h those xed in the original mo del, (3 ) . (Note that, in [ 37 ℄, a dieren t limit w as emplo y ed, namely , n → ∞ with lim n →∞ n ( c ) /n ( v ) = C 1 ∈ (0 , ∞ ) and lim n →∞ n ( f ) /n = 0 ; this is appropriate for exat asymptotis, but not for sim ulation, b eause the asymptoti tilting fator to b e used in the latter then do es not feel the foreign an tigens.) Lemma 2 L et f b e the density of Exp (1 /τ ) (i.e., f ( τ ) = e − τ / τ / τ ), and ψ ( t ) := E ( e tW ) = Z ∞ 0 exp tw ( τ ) f ( τ ) dτ = 1 ¯ τ Z ∞ 0 exp t exp( − 1 /τ ) τ − τ ¯ τ dτ (31) b e the moment-gener ating funtion of W 1 . Under the assumptions of Se t. 4.3 , the unique solution ϑ ∗ of g act m = m ( c ) m q z ( c ) d d t log ψ ( t ) t = qz ( c ) ϑ + m ( v ) m q z ( v ) d d t log ψ ( t ) t = qz ( v ) ϑ + 1 m z ( f ) d d t log ψ ( t ) t = z ( f ) ϑ (32) is the unique asymptoti al ly eient tilting p ar ameter for LD simulation of P n ( A ) . Pr o of Clearly , P n satises the assumptions of Set. 4.1 . Note, in partiular, that ψ ( t ) < ∞ for all t ∈ R sine W is b ounded ab o v e and b elo w, and so Λ ( ϑ ) = lim n →∞ log E ( e ϑG n ( z ( f ) ) /n ) = m ( c ) m log ψ ( q z ( c ) ϑ ) + m ( v ) m log ψ ( q z ( v ) ϑ ) + 1 m log ψ ( z ( f ) ϑ ) < ∞ (33) for all ϑ ; hene, the Gärtner-Ellis theorem holds b y Lemma 1 . T o v erify the remaining assumptions of Prop. 1, reall from Se. 4.1 that Λ ( ϑ ) is dieren tiable (with on tin uous deriv ativ e) on all of R . The b ounds on g act /m lead to Λ ′ (0) = E G ( z ( f ) m < g act m < M w (1) m = lim ϑ →∞ Λ ′ ( ϑ ) . (34) Λ is stritly on v ex (sine (d 2 / d t 2 ) log ψ ( t ) is the v ariane of W t , the tilted v ersion of W (f. [ 2, Prop. XI I.1.1℄), whi h is p ositiv e sine W and hene W t is nondegenerate). Eq. (34 ) th us en tails that Λ ′ ( ϑ ) = g act /m has a unique solution ϑ ∗ , whi h is p ositiv e (and learly satises (V2)). As a 14 onsequene, g act /m is a dominating p oin t of A , whi h is a rare ev en t sine 0 < I ( g act /m ) < ∞ (b y Λ (0) = 0 together with (34 ) and (12 ) ; f. Fig. 4, left). Finally , A is a on tin uit y set of b oth I and I + h ϑ ∗ , . i simply b eause I and h ϑ ∗ , . i are on tin uous at g act /m , and A = A ◦ . Realising that the righ t-hand side of (32 ) equals Λ ′ ( ϑ ) (see also Eq. (20) in [37 ℄), one obtains the laim from Prop. 1. ⊓ ⊔ 0 0 . 01 0 . 02 0 . 03 0 . 04 0 . 05 0 . 06 0 0 . 03 0 . 06 ϑ Λ ( ϑ ) aϑ I ( a ) ϑ ⋆ 0 0 . 004 0 . 008 0 . 012 0 . 016 300 600 900 I ( g act /m ) g act Fig. 4 The um ulan t-generating funtion Λ (left) and the rate funtion I (righ t) for the T-ell mo del (3). The slop e of the straigh t line in the left panel is a = g act /m , where g act = 800 and m = 1551 . A t ϑ ∗ , aϑ − Λ ( ϑ ) assumes its maxim um, I ( a ) (f. (10 ) (12 )). The solution of (32 ) is readily alulated n umerially . The funtion Λ , and the resulting rate funtion I , are sho wn in Fig. 4. As desrib ed in Set. 4.2 , w e no w tilt the densit y f of the T j with ϑ ∗ aording to Eq. (25 ) . This yields three dieren t densities ˜ f ϑ ∗ α , dep ending on the w eigh ting fators α ∈ { q z ( c ) , q z ( v ) , z ( f ) } , namely ˜ f ϑ ∗ α ( τ ) = exp( αϑ ∗ w ( τ )) f ( τ ) ψ ( αϑ ∗ ) = 1 ¯ τ exp αϑ ∗ exp( − 1 /τ ) τ − τ ¯ τ ψ ( αϑ ∗ ) . (35) As disussed in Set. 4.2, this is not the densit y of an y kno wn standard distribution (let alone an exp onen tial one), and sim ulating from it requires n umerial in tegration (whi h is w ell-b eha v ed sine the ˜ f ϑ ∗ α are n umerially w ell-b eha v ed), and disretisation and tabulation of the resulting distribution funtions ˜ F ϑ ∗ α , follo w ed b y lo oking up the solution ˜ τ ϑ ∗ of ˜ F ϑ ∗ α ( ˜ T ϑ ∗ ) = U for U ∼ Uni [0 , 1] , to nally yield αW ϑ ∗ via αW ϑ ∗ = αw ( ˜ T ϑ ∗ ) . Sear hing the table w ould b e the sp eed- (or preision-) limiting step, requiring O (lo g D ) op er- ations if D is the n um b er of disretisation steps. This an b e remedied b y applying the so-alled alias metho d to qui kly generate random v ariables aording to the disretised probabilit y distri- bution. F or a desription of the metho d, w e refer the reader to [19, pp. 2527℄, [16 ℄, or [24 , p. 248℄. Let us just summarise here that, after a prepro essing step, whi h is done one for a giv en distri- bution, the metho d only requires one Uni [0 , 1] random v ariable together with one m ultipliation, one uto and one subtration (or t w o Uni [0 , 1] random v ariables together with one m ultipliation, one uto and one omparison, dep ending on the implemen tation) to generate one realisation of ˜ T ϑ ∗ , regardless of D (in partiular, it do es without sear hing altogether). W e no w ha v e ev erything at hand to form ulate the algorithm to sim ulate (realisations of ) G ( z ( f ) ) of (3). (F or notational on v eniene, w e will not distinguish b et w een random v ariables and their realisations here). Algorithm 1 ompute ϑ ∗ by solving Eq. (32 ) numeri al ly 15 alulate the tilte d densities ˜ f ϑ ∗ α , α ∈ { q z ( c ) , q z ( v ) , z ( f ) } , via (35 ) for i=1 til l sample size N do for every summand j of (3) gener ate a sample ( ˜ T ϑ ∗ j ) ( i ) a or ding to its density ˜ f ϑ ∗ α ( j ) with the help of the alias metho d (her e, the upp er index ( i ) is adde d to r ee t sample i , and α ( j ) is the weighting fator of the sum to whih j b elongs) alulate G ( z ( f ) ) ( i ) = m ( c ) X j =1 q z ( c ) w ( ˜ T ϑ ∗ j ) ( i ) ! + m ( c ) + m ( v ) X j = m ( c ) +1 q z ( v ) w ( ˜ T ϑ ∗ j ) ( i ) ! + z ( f ) w ( ˜ T ϑ ∗ m ( c ) + m ( v ) +1 ) ( i ) alulate the indi ator funtion times the r eweighting fator (i.e., the i -th summand in Eq. (15 ) ) if ( G ( z ( f ) )) ( i ) ≥ g act then R ( i ) = m Y j =1 f α ( j ) (( ˜ T ϑ ∗ j ) ( i ) ) ˜ f ϑ ∗ α ( j ) (( ˜ T ϑ ∗ j ) ( i ) ) else R ( i ) = 0 end if end for alulate \ P ϑ ∗ P m ( A ) N = P N i =1 R ( i ) N , as estimate of P ( G ( z ( f ) ) > g act ) . 4.4 Extension to v ariable op y n um b ers Let us no w onsider the extended mo del (4), in whi h the op y n um b ers are themselv es random v ariables. This is also o v ered b y the large deviation theory presen ted ab o v e; in partiular, Lemma 1 again applies if the Y ( k ) ℓ in (23) are iden tied with Z ( c ) j W j or Z ( v ) j W j , resp etiv ely . The global tilting fator ϑ ∗ is, in the usual w a y , alulated as the solution of Λ ′ ( ϑ ) = g act /m , where Λ ( ϑ ) is as in (33) with ψ ( q z ( k ) ϑ ) = E ( e qz ( k ) ϑW ) replaed b y E ( ψ ( q Z ( k ) ϑ )) = E ( e qZ ( k ) ϑW ) , k ∈ { c, v } ; see Eq. (20) in [37 ℄. Ho w ev er, the ob jet of tilting no w is the joint distribution of W and Z ( c ) (or Z ( v ) , resp etiv ely), that is, d F ( τ )d H ( k ) ( z ) reeiv es the rew eigh ting fator exp( q ϑz w ( τ )) , where F and H ( k ) denote the measures of T and Z ( k ) , k ∈ { c , v } , resp etiv ely . This in tro dues dep endenies b et w een op y n um b ers and stim ulation rates. The resulting bivariate sim ulation task is ostly and ma y oset some of the eieny gain obtained b y tilting. If, ho w ev er, the Z ( k ) are losely p eak ed around their means (as is the ase for our hoie of parameters), the follo wing h ybrid pro edure turns out to b e b oth pratial and fast: Dra w the Z ( k ) from their original (un tilted, binomial) distributions; and sim ulate a tilted v ersion of q W , denoted b y ( q W ) ϑ ∗ , b y rew eigh ting the original densit y of q W with exp( q ϑ ∗ E ( Z ( k ) ) W ) , irresp etiv e of the atual v alue of Z . Clearly , this metho d is not asymptotially eien t, but it is a v alid imp ortane sampling metho d that turns out to ompare w ell with the ideal pro edure used for the xed op y n um b ers (see Se. 5.1.3 ). 5 Results Let us no w presen t the results of our sim ulations in t w o steps. W e rst in v estigate the p erformane of the metho d, and then use it to gain more insigh t in to the underlying phenomenon of statistial reognition. 16 5.1 P erformane of the sim ulation metho d W e will examine the p erformane of the imp ortane-sampling metho d in three resp ets: w e will ompare it to simple sampling (the previously-used sim ulation metho d) and to the results of exat asymptotis (the previously-used analyti metho d); nally , w e will quan tify the eieny in terms of the relativ e error (and th us return to the theory of Set. 3.2 ). In an y ase, w e will onsider P ( G ( z ( f ) ) ≥ g act ) as a funtion of g act (and for v arious v alues of the parameter z ( f ) ). Of ourse, this probabilit y is just one min us the distribution funtion of G ( z ( f ) ) ; in imm unobiology , the orresp onding graph is kno wn as the ativ ation urv e. Ev aluating this graph b y LD sim ulation requires, for ea h v alue of g act to b e onsidered, a fresh sample, sim ulated with its individual tilting fator ϑ ∗ (reall that this dep ends on g act via (32 ) ). A t rst sigh t, this lo oks lik e an enormous disadv an tage relativ e to simple sampling, where no threshold needs to b e sp eied in adv ane; rather, the outomes of the sim ulation diretly yield an estimate o v er the en tire range of the ativ ation urv e. Ho w ev er, it will turn out that this disadv an tage is oset man y times b y the sp ei eieny of hitting the rare ev en ts in LD sampling. (There is ro om for impro v emen t: the samples that do not hit a giv en rare ev en t ould b e used to impro v e the estimates of the more lik ely ev en ts.) 5.1.1 Comp arison with simple sampling Clearly , b oth the simple-sampling and the imp ortane-sampling estimates are un biased and on- v erge to the true v alues as N → ∞ . It is therefore no surprise that they yield pratially iden tial results wherev er they an b e ompared and this yields a rst qui k onsisteny he k for our metho d. This is demonstrated in Fig. 5, whi h sho ws simple sampling (SS) and imp ortane sampling (IS) ativ ation urv es, ea h for z ( f ) = 10 00 and z ( f ) = 20 00 . F or SS, N = 1 . 3 ∗ 10 8 samples, ( G ( z ( f ) )) ( i ) , 1 ≤ i ≤ N , w ere generated altogether for ev ery graph, whereas for IS, N = 1000 0 samples w ere generated for ev ery threshold v alue onsidered (from g act = 100 to g act = 1000 in steps of 50 ), i.e. 1 . 9 ∗ 10 5 samples altogether. Bey ond g act = 450 and g act = 800 (for z ( f ) = 1000 and z ( f ) = 2000 , resp etiv ely), no estimates ould b e obtained via SS due to the lo w probabil- ities in v olv ed, whereas with IS, it is easy to get b ey ond g act = 900 in either ase, although the probabilities an get do wn to 10 − 20 (note, ho w ev er, that this far end of the distribution is no longer biologially relev an t). In terms of run time, determining an ativ ation urv e (o v er its en tire range) b y SS to ok 48 hours of CPU time (In tel P en tium M 1.4 GHz 512MB RAM), whereas IS required only ab out 2 min utes (in the threshold regime where the metho ds are omparable), that is, a sp eedup b y a fator of nearly 1500 is a hiev ed. W e also applied our metho d to the extended mo del (4) with binomially distributed op y n um- b ers. Figure 6 sho ws the sim ulation results for t w o v alues of z ( f ) , ea h for SS and IS. Again, the urv es agree, as they m ust. As to run time, it to ok ab out 130 hours to generate the 2 ∗ 1 0 7 samples for SS, whereas for IS it to ok 10 min. to generate the 9 . 5 ∗ 10 4 samples. 5.1.2 Comp arison with exat asymptotis A pillar of the previous analysis of Zin t et al. [ 37 ℄ (and its preursor BRB [33 ℄) has b een so-alled exat asymptotis. This is a renemen t of large deviation theory whi h yields estimates for the probabilities P n ( A ) themselv es, rather than just their exp onen tial dea y rates obtained via the LDP in Def. 1. With standard large deviation theory (and our sim ulation metho d), it shares the tilting parameter whi h is alulated aording to Eq. (32 ); for more details, w e refer to [37 ℄. A omparison of IS sim ulation with exat asymptotis is also inluded in Fig. 5 . F or small v alues of g act , exat asymptotis is sligh tly impreise. This is due to the asymptoti nature ( n → ∞ ) of the metho d, whi h yields more preise results in the v ery tail of the distribution, where the deviations are truly large. Note that, although our tilting fators agree with those in exat asymptotis, rare ev en t sim ulation do es not suer from this auray problem sine, due to the rew eigh ting, it is 17 10 − 20 10 − 12 10 − 06 0.01 100 300 600 900 P ( G ( z ( f ) ) ≥ g act ) g act IS for z ( f ) = 0 SS for z ( f ) = 1000 IS for z ( f ) = 1000 LDT for z ( f ) = 1000 SS for z ( f ) = 2000 IS for z ( f ) = 2000 LDT for z ( f ) = 2000 Fig. 5 Estimates of the ativ ation urv e, P ( G ( z ( f ) ) ≥ g act ) , in the basi mo del (3 ) for z ( f ) = 1000 and z ( f ) = 2000 , as w ell as for the self ba kground ( z ( f ) = 0 ), on logarithmi sale. The probabilities w ere estimated indep enden tly with simple sampling (SS), imp ortane sampling (IS), and exat asymptotis based on large deviation theory (LDT) as used in [37 ℄. F or IS, 19 v alues of g act w ere onsidered (from 100 to 1000 in steps of 50 ), and N = 10000 samples w ere generated for ea h v alue (i.e., 1 . 9 ∗ 10 5 samples altogether), whereas for the SS sim ulation, N = 1 . 3 ∗ 10 8 samples w ere used o v er the en tire range. The SS urv es end at g act = 400 and g act = 800 , resp etiv ely , b eause larger v alues w ere not hit in the giv en sample. The IS and SS graphs agree p erfetly un til the SS sim ulation la ks preision. F or larger threshold v alues, w e see a p erfet agreemen t of the IS and LDT graphs. Note the general feature that, for threshold v alues that are not to o small, the ativ ation probabilit y in the presene of foreign an tigens is sev eral orders of magnitude larger than the self ba kground, i.e. Eq. (6) is satised. alw a ys a v alid imp ortane sampling s heme that yields un biased estimates for ev ery nite n ; the nite-size eets will only manifest themselv es as a ertain loss of eieny , as will b e seen b elo w. 5.1.3 Asymptoti eieny and r elative err or In order to in v estigate the relativ e error of ( \ P P ϑ ∗ n ( A )) N , w e rst note that the v ariane of the estimator is giv en b y V \ P P ϑ ∗ n ( A ) N = 1 N V \ P P ϑ ∗ n ( A ) 1 = 1 N E h 1 { ( T ϑ ∗ n ) (1) ∈ A } d P d P ϑ ∗ n ( T ϑ ∗ n ) (1) − P n ( A ) 2 i , (36) where w e ha v e used (15 ) for N = 1 . V ( \ P P ϑ ∗ n ( A )) 1 an b e estimated via the giv en n um b er N of samples in a single sim ulation run, i.e., as the sample v ariane b V \ P P ϑ ∗ n ( A ) 1 = 1 N − 1 N X i =1 1 { ( t ϑ ∗ n ) ( i ) ∈ A } d P d P ϑ ∗ n ( t ϑ ∗ n ) ( i ) − \ P P ϑ ∗ n ( A ) N 2 , (37) 18 10 − 14 10 − 10 10 − 06 0 . 01 300 600 900 P ( G ( z ( f ) ) ≥ g act ) g act SS for z ( f ) = 2500 IS for z ( f ) = 2500 SS for z ( f ) = 1500 IS for z ( f ) = 1500 Fig. 6 Sim ulation of P ( G ( z ( f ) ) ≥ g act ) in the extended mo del (4), for z ( f ) = 1500 and z ( f ) = 2500 . The probabilities w ere estimated indep enden tly with simple sampling, and with imp ortane sampling at 19 dieren t threshold v alues (from 100 to 1000 in steps of 50 ). F or IS, 9 . 5 ∗ 10 4 samples w ere generated ( 5000 p er threshold); for SS, 2 ∗ 10 7 samples w ere used. No estimates are obtained with SS at thresholds b ey ond 600 or 920, resp etiv ely , in analogy with the situation in Fig. 5 . where the ( t ϑ ∗ n ) ( i ) are no w onsidered as realisations of ( T ϑ ∗ n ) (1) . W e an th us estimate the squared relativ e error as c η 2 N ( P ϑ ∗ n , A ) = 1 N b V \ P P ϑ ∗ n ( A ) 1 \ P P ϑ ∗ n ( A ) N 2 . (38) F or simple sampling, one pro eeds in the ob vious analogous w a y (without tilting and rew eigh ting). In line with the limit disussed in Se. 4.3 , w e no w onsider G n ( z ( f ) ) for system sizes n = n i , where n i = n ( c ) i + n ( v ) i + n ( f ) i , 0 ≤ i ≤ 10 , and w e ho ose n ( α ) i = im ( α ) , α ∈ { c, v, f } , for 1 ≤ i ≤ 10 , as w ell as n ( c ) 0 = m ( c ) / 2 , n ( v ) 0 = m ( v ) / 2 , and n ( f ) 0 = m ( f ) (i.e., w e simply `m ultiply' the system, exept for i = 0 , whi h orresp onds to `half ' a system exept for the foreign p eptide, whi h annot b e split in to t w o). W e then sim ulate P ( G n i ( z ( f ) ) ≥ g act n i /m ) for t w o v alues of z ( f ) and a xed v alue of g act with our imp ortane sampling metho d, as sho wn in Fig. 7 . Ob viously , the (estimated) probabilities dea y to zero at an exp onen tial rate with inreasing n , as they m ust b y their LDP . In on trast, the (estimated) squared RE only inreases linearly this ev en go es b ey ond the predition of the theory (asymptoti eieny only guaran tees a sub exp onen tial inrease). So far, w e ha v e onsidered the n -dep endene of the metho d for a xed v alue of g act , in the ligh t of the a v ailable asymptoti theory . F or the pratial sim ulation of the giv en T-ell problem, w e no w tak e the giv en system size n = m and n umerially in v estigate the relativ e error as a funtion of g act . Here, the exp onen tial dea y of P ( G ( z ( f ) ) ≥ g act ) as a funtion of g act is deisiv e, whi h w e ha v e already observ ed in Fig. 5 , and whi h go es together with the at-least-linear inr e ase of I with g act (reall that I is on v ex, and see Fig. 4). Fig. 8 sho ws the relativ e error of b oth SS and IS. It do es not ome as a surprise that, again, IS do es extremely w ell and b eats the exp onen tial 19 10 − 50 10 − 40 10 − 30 10 − 20 1 − 10 1 0 2 4 6 8 10 P ( G n ( z ( f ) ) ≥ 400 n m ) n/m z ( f ) = 1000 z ( f ) = 2000 0 . 0005 0 . 001 0 . 0015 0 . 002 0 . 0025 0 . 003 0 2 4 6 8 10 c η 2 N n/m z ( f ) = 1000 z ( f ) = 2000 Fig. 7 Imp ortane sampling sim ulations for P ( G n ( z ( f ) ) ≥ g act n/m ) for n = n i , 0 ≤ i ≤ 10 , for g act = 400 and t w o v alues of z ( f ) . Left: estimate of the probabilit y (note that the v ertial axis is on logarithmi sale). Righ t: estimated squared RE. dea y of the probabilities: whereas, on the log sale of the v ertial axis, the squared RE of SS gro ws roughly linearly , it remains more or less onstan t for IS. (The v ery lo w squared RE of the simple sampling graphs for lo w thresholds in the righ t panel is due to the fat that the probabilit y to rea h this threshold is quite high and the h uge sample of N = 1 . 3 · 1 0 8 on tributes to estimating it, that is, the sample sizes are not omparable. A simple sampling sim ulation run with the total sample size of a orresp onding IS sim ulation (i.e., N = 1 0000 times the n um b er of steps on tained in the in terv al onsidered) results in higher relativ e errors than for imp ortane sampling ev en for the lo w threshold v alues (left panel). W e w ould lik e to note, ho w ev er, that the run time of simple sampling for these small sample sizes is shorter than the run time for IS, ev en if one do es not oun t the o v erhead required to get the tilting parameters for imp ortane sampling.) 0 . 001 0 . 01 1 100 300 600 900 c η 2 N g act IS-RE for z ( f ) = 1000 IS-RE for z ( f ) = 2000 SS-RE for z ( f ) = 1000 SS-RE for z ( f ) = 2000 10 − 7 0 . 0001 0 . 01 100 300 600 900 c η 2 N g act IS-RE for z ( f ) = 1000 IS-RE for z ( f ) = 2000 SS-RE for z ( f ) = 1000 SS-RE for z ( f ) = 2000 Fig. 8 Estimated squared RE for simple sampling ( N = 10000 times the n um b er of steps on tained in the onsidered in terv al (left), N = 1 . 3 ∗ 10 8 (righ t)), and imp ortane sampling ( N = 10000 p er threshold v alue in either panel) sim ulations of P ( G ( z ( f ) ) > g act ) of the basi mo del, Eq. (3). Note that the v ertial axis is on logarithmi sale. Figure 9 sheds more ligh t on the b eha viour of the relativ e error of the IS sim ulation. It sho ws the squared RE for 6 distint z ( f ) -v alues and rev eals the nite-size eets. The w a v e-lik e b eha viour for larger z ( f ) is due to the fat that, for v ery lo w threshold v alues, there is no real need for tilting, b eause the original distribution P n is already lose to optimal and the tilting fator is v ery small. F or inreasing thresholds, substan tial tilting is required, but there are still visible deviations from 20 0 0.002 0.004 0.006 100 250 500 750 1000 c η 2 N g act z ( f ) = 250 z ( f ) = 500 z ( f ) = 1000 z ( f ) = 1500 z ( f ) = 2000 z ( f ) = 2500 Fig. 9 Estimated squared RE of our IS estimate, for v arious frequenies z ( f ) of the foreign an tigen. Details are as in Fig. 8, but no w the v ertial axis is on linear sale. the n → ∞ limit (as already disussed in the on text of Fig. 5), so the tilted distributions are not optimal. This pro dues the h ump in the squared RE urv es, whi h is more pronouned for larger z ( f ) v alues b eause, for the ase n = m onsidered here, the foreign an tigens ome as a single term that ma y stand out. F or large g act , nally , one gets lose enough to the limit, and the exp eted sub-exp onen tial inrease sets in (in our ase, it is, in fat, roughly linear). Nev ertheless, it should b e lear that, in spite of the sligh t non-optimalit y at small threshold v alues, our tilted distributions still yield a far lo w er squared RE than do es simple sampling. A v ery similar piture emerges for the extended mo del; surprisingly , the relativ e error is no larger than in the basi mo del, although the ad ho sim ulation metho d used here is not asymptotially eien t (see Se. 4.4 ; data not sho wn). 5.2 Analysis of the T-ell mo del In this Setion, w e use our sim ulation metho d to obtain more detailed insigh t in to the phenomenon of statistial reognition in the T-ell mo del. As disussed b efore, the task is to disriminate one foreign an tigen t yp e against a noisy ba kground of a large n um b er of self an tigens. W e already kno w from Fig. 5 that, for threshold v alues that are not to o small, the ativ ation probabilit y in the presene of foreign an tigens is sev eral orders of magnitude larger than the ativ ation probabilit y of the self-ba kground, i.e. Eq. (6) is satised. As disussed in [37 ℄, this distintion relies on z ( f ) > z ( c ) , z ( v ) what happ ens is that larger op y n um b ers of the foreign an tigen thi k en the tail of the distribution of G ( z ( f ) ) (without hanging its mean), so that the threshold is more easily surpassed. The self-nonself distintion ma y , aording to this mo del, b e roughly desrib ed as follo ws. F or a giv en an tigen (foreign or self ), nding a highly-stim ulating T-ell reeptor is a rare ev en t; but if it o urs to a foreign an tigen, it o urs man y times sim ultaneously sine there are n umerous opies, whi h all on tribute the same large signal, sine all reeptors of the T-ell in v olv ed are iden tial; the resulting stim ulation rate is th us high. In on trast, if it is a self an tigen that nds a highly-stim ulating reeptor, the eet is less pronouned due to the smaller op y n um b ers. In this sense, the to y mo del explains the distintion solely on the basis of op y n um b ers; but see the Disussion for more sophistiated eets that alleviate this requiremen t. F ollo wing these in tuitiv e argumen ts, w e no w aim at a more detailed piture of ho w the self ba kground lo oks, and ho w the foreign t yp e stands out against it. T o in v estigate this, it is useful 21 rate \ g act 0 100 250 500 1000 v ariable 66 . 6 74 . 9 77 . 1 78 . 8 80 . 0 constitutive 22 . 2 59 . 2 277 . 7 590 . 6 1160 rate \ g act 0 100 250 500 1000 v ariable 12 . 7 13 . 9 14 . 5 14 . 9 15 . 1 constitutive 23 . 1 35 . 6 88 . 8 134 . 9 191 . 3 T able 1 Sample means (left) and sample standard deviations (righ t) of the histograms in Fig. 10 (left) and Fig. 11 (i.e., the self-only ase). rate \ g act 0 100 250 500 1000 v ariable 65 . 9 74 . 1 74 . 2 76 . 2 78 . 4 constitutive 21 . 8 55 . 9 129 . 5 270 . 4 821 . 1 foreign 0 . 9 4 . 0 184 . 8 279 . 6 302 . 2 rate \ g act 0 100 250 500 1000 v ariable 12 . 7 14 . 1 13 . 9 14 . 2 14 . 7 constitutive 22 . 4 42 . 0 90 . 4 109 . 1 163 . 7 foreign 6 . 7 18 . 5 112 . 2 54 . 5 39 . 2 T able 2 Sample means (left) and sample standard deviations (righ t) of the histograms in Fig. 10 (righ t) and Fig. 12 (i.e., the ase with foreign an tigens). to onsider the histograms of the total onstitutiv e, v ariable, and foreign stim ulation rates, i.e., the on tributions of the onstitutiv e sum, the v ariable sum, and the individual foreign term in the sum (3), either for all samples or for the subset of samples for whi h G ( z f ) ≥ g act , for v arious g act . Sine this requires a higher resolution (and th us larger sample size) than the alulation of the ativ ation probabilities alone, su h analysis w ould b e pratially imp ossible with simple sampling. With IS, w e again generated 10000 samples p er g act v alue, from whi h b et w een 30 and 70 p eren t turned out to rea h the threshold. 0 1000 2000 3000 0 100 200 300 stim ulation rate onstitutiv e v ariable 0 3000 6000 9000 0 100 200 300 stim ulation rate onstitutiv e foreign v ariable Fig. 10 Histograms of the total stim ulation rates of v ariable, onstitutiv e, and foreign an tigens, for z ( f ) = 0 (left) and z ( f ) = 1000 (righ t), in the basi mo del (3), when all samples are inluded. Sample size is 10000, and the v ertial axis holds the n um b er of samples whose total onstitutiv e (v ariable, foreign) stim ulation rates fall in to giv en in terv als. Note that the saling of the v ertial axis v aries aross diagrams. Figure 10 sho ws the resulting histograms when all samples are inluded, and Figs. 11 and 12 sho w the histograms for the subset of samples that ha v e surpassed four represen tativ e threshold v alues, without and with foreign an tigen. T ables 1 and 2 summarise these results in terms of means and standard deviations. Finally , Fig. 13 sho ws the orresp onding t w o-dimensional statistis for all pairs of v ariable, onstitutiv e, and foreign stim ulation rates, again for v arious threshold v alues. (Figs. 11 13 are based on the outome of imp ortane sampling without r eweighting ; normalising b y the n um b er of "suessful" samples w ould result in an estimate of the onditional distribution, b eause the rew eigh ting fators anel out.) Let us start with the situation without foreign an tigens, as displa y ed in Figs. 10 (left) and 11 as w ell as T able 1. This already illustrates the fundamen tal dierene b et w een v ariable and onstitutiv e an tigens. Judging from the large n um b er ( m ( v ) = 15 00 ) of individual terms in the sum at lo w op y n um b er ( z ( v ) = 50 ), the v ariable stim ulation rate is exp eted to b e appro ximately 22 0 250 500 0 50 100 150 200 250 stim ulation rate onstitutiv e v ariable 0 200 400 600 0 200 400 600 800 stim ulation rate onstitutiv e v ariable 0 200 400 600 0 300 600 900 1200 1500 stim ulation rate onstitutiv e v ariable 0 200 400 0 400 800 1200 1600 2000 stim ulation rate onstitutiv e v ariable Fig. 11 Histograms of the total stim ulation rates of v ariable and onstitutiv e an tigens, for z ( f ) = 0 , in the basi mo del (3 ) , for samples that rea h a giv en threshold v alue ( g act = 10 0 (upp er left), g act = 25 0 (upp er righ t), g act = 500 (lo w er left), g act = 1000 (lo w er righ t)). Sample size is 10000, and the v ertial axis holds the n um b er of samples that rea h g act and whose total onstitutiv e (v ariable, foreign) stim ulation rates falls in to giv en in terv als. Note that the saling of b oth axes v aries aross diagrams. normally distributed and fairly losely p eak ed around its mean at least as long as no restrition on G ( z ( f ) ) is in v olv ed and, as the Figure sho ws, this feature p ersists when G ( z ( f ) ) > g act , pratially indep enden tly of the threshold in v olv ed. So, the v ariable an tigens form a kind of ba kground that p oses no diult y to foreign-self distintion: it is not v ery noisy , and it do es not hange with the threshold. In on trast, the distribution of the onstitutiv e ativ ation rates is wider; this is due to the large op y n um b ers ( z ( c ) = 500 ), the eet of whi h is not omp ensated b y the smaller n um b er of terms, m ( c ) = 50 . F urthermore, the normal appro ximation is not exp eted to b e partiularly go o d for the onstitutiv e an tigens giv en the extreme asymmetry of the W -distribution (see Fig. 3), the en tral limit theorem will not a v erage out the deviations at only m ( c ) = 50 . In partiular, the distribution remains asymmetri. With inreasing threshold, this distribution mo v es to the righ t. The reason for this is that, in order to rea h an inreasing g act , the tail ev en ts of the onstitutiv e or the v ariable sum or b oth m ust b e used, but it is easier (that is, more probable) to use the onstitutiv e one b eause it on tains more at ypial ev en ts. In the language of large deviation theory , this is an example of the general priniple that large deviations are alw a ys done in the the least unlik ely of all the unlik ely w a ys [8 , Ch. I℄. In the language of biology , the onstitutiv e an tigens are the problem of foreign-self distintion: due to their high op y n um b ers and inomplete a v eraging, utuations p ersist that o asionally indue an imm une resp onse ev en in the absene of foreign an tigens. This o urs if a T-ell reeptor happ ens to t partiularly w ell to one, or a n um b er of, onstitutiv e an tigen t yp es on an APC; due to their large op y n um b ers, 23 0 500 1000 1500 2000 0 100 200 300 stim ulation rate onstitutiv e foreign v ariable 0 200 400 0 200 400 stim ulation rate onstitutiv e foreign v ariable 0 200 400 600 0 300 600 900 stim ulation rate onstitutiv e foreign v ariable 0 200 400 0 400 80 0 1200 1600 stim ulation rate onstitutiv e foreign v ariable Fig. 12 Histograms of the total onstitutiv e, v ariable and foreign stim ulation rates for z ( f ) = 1000 in the basi mo del (3). Sample size is 10000, and the v ertial axis holds the n um b er of samples that rea h the threshold g act and whose total onstitutiv e (v ariable, foreign) stim ulation rate falls in to a giv en in terv al, for g act = 100 (upp er left), g act = 250 (upp er righ t), g act = 500 (lo w er left), g act = 1000 (lo w er righ t). The maximal stim ulation rate for the foreign an tigens is z ( f ) w ( 1) = 367 . 9 . Note that the saling of b oth axes v aries aross diagrams. these few highly-stim ulating t yp es are then suien t to surpass the threshold (in on trast, sev eral highly-stim ulating t yp es w ould b e required for the v ariable an tigens to eliit a reation, whi h is to o improbable). Let us no w turn to the piture with foreign an tigen presen t (Figs. 10 (righ t), 12 , 13 , and T able 2 ). One salien t feature here is that the v ariable stim ulation rate b eha v es exatly as in the self-only ase: losely p eak ed around a small mean, un hanged when { G ( z ( f ) ) > g act } is imp osed. The piture is th us dominated b y the in terpla y of onstitutiv e and foreign t yp es. In line with Fig. 5, the situation is similar in the ase without restrition on G ( z ( f ) ) (Fig. 10 , righ t) and the ase when G ( z ( f ) ) ≥ 100 (Fig. 12 , upp er left). In partiular, the foreign stim ulation rate is losely p eak ed at 0 ; only the onstitutiv e ba kground has mo v ed sligh tly to the righ t, exatly as in the self-only ase. F or g act = 250 (Fig. 12 , upp er righ t), where, aording to Fig. 5 , foreign-self distintion sets in, the foreign stim ulation rate b eomes prominen t: the righ t bran h of the W -distribution no w b eomes p opulated, and the asso iated stim ulation rates are large due to the large op y n um b ers z ( f ) in v olv ed. Nev ertheless, for g act = 250 , the foreign stim ulation rate is lose to 0 in a sizable fration of the ases in whi h an imm une reation o urs here, the reation is brough t ab out b y the onstitutiv e ba kground, whi h mo v es to the righ t just as in the self-only ase (but less pronouned). Fig. 13 sho ws that the onstitutiv e and foreign stim ulation rates are, indeed, negativ ely orrelated: as is to b e exp eted, lo w foreign rates are omp ensated b y high onstitutiv e rates and vie v ersa 24 (in on trast, the v ariable ba kground hardly orrelates with either the onstitutiv e or the foreign stim ulation rate). As in the self-only ase, therefore, the lev el of un w an ted ativ ation (self-only or mainly self, without appreiable foreign ativ ation) is set b y the tail b eha viour of the onstitutiv e ba kground. Ho w ev er, if g act is inreased further (Fig. 12 , lo w er left), ev ery T ell b ey ond the threshold displa ys high stim uli for the foreign an tigen, their distribution shifting ev en further to the righ t and onen trating near the maximal stim ulation rate giv en b y the maxim um of the funtion w of Eq. ( 1), more preisely , b y z ( f ) w (1) . This maxim um an, of ourse, not hange b y imp osing restritions on G ( z ( f ) ) ; th us, an y further inrease of g act (Fig. 12 , lo w er righ t) m ust then b e mat hed b y the b y no w familiar shift of the onstitutiv e ba kground. (This last panel is, ho w ev er, less biologially realisti sine the probabilities in v olv ed are to o small to b e relev an t after all, with ab out 10 7 dieren t T-ell t yp es, threshold v alues that yield ativ ation probabilities far b elo w 10 − 7 ev en in the presene of foreign an tigens oer no imm une protetion.) A further illustration of the onset of self-nonself distintion is presen ted in Fig. 14 . Here w e onsider P G ( z ( f ) ) − z ( f ) W n ( c ) + n ( v ) +1 > g act | G ( z ( f ) ) > g act = P G ( z ( f ) ) − z ( f ) W n ( c ) + n ( v ) +1 > g act P G ( z ( f ) ) > g act , (39) i.e., the probabilit y that, in a T-ell that is ativ ated in the presene of foreign an tigen, the self omp onen t alone w ould ha v e b een suien t for the ativ ation. F rom z ( f ) = 1000 on w ards, this probabilit y dereases to 0 qui kly with inreasing g act . Put dieren tly , in large parameter regions, the foreign an tigens do indeed mak e the dierene, whi h is the deisiv e feature of self-nonself distintion. 6 Conlusion and outlo ok W e ha v e established here a metho d of LD sampling that allo ws the on v enien t sim ulation of the rare ev en ts relev an t to statistial reognition in the imm une system. Th us a more thorough in v estigation of these ev en ts ould b e arried out. But this is only a rst step, and the goal for future w ork is to use this or related metho ds to in v estigate biologially realisti mo dels. Indeed, the to y mo del onsidered here, whi h relies solely on distintion b y op y n um b ers, do es serv e the aim to illustrate that distintion against a noisy ba kground is, at all, p ossible, ev en without an in trinsi dierene b et w een self and nonself, and ho w this is related to the rare ev en ts in the tail of the ba kground distribution. Ho w ev er, biologially realisti mo dels ha v e to tak e in to aoun t tolerisation me hanisms that mak e the T- ells less resp onsiv e to self an tigens. One imp ortan t su h me hanism is so-alled ne gative sele tion . Negativ e seletion o urs during the maturation phase of y oung T-ells in the th ym us, b efore they are released in to the b o dy . In a pro ess similar to the one desrib ed b y the to y mo del, they are onfron ted with APCs that presen t mixtures of v arious self an tigens, and those T-ells whose ativ ation rate surpasses a th ymi ativ ation threshold g thy < g act are eliminated. When they are later, after lea ving the th ym us, onfron ted with mixtures of self and foreign an tigens, the stim ulation rates emerging from self and foreign are no longer i.i.d. (the self ones are biased to w ards smaller v alues and p ossibly negativ ely orrelated). In fat, a simple mo del for negativ e seletion w as already desrib ed in BRB [33 ℄, and sho wn to drastially redue the self ba kground, so that foreign an tigens do no longer require elev ated op y n um b ers to b e deteted. More sophistiated mo dels of negativ e seletion ha v e b een form ulated e.g. in [32 ℄. Ho w ev er, their sim ulation still a w aits the dev elopmen t of adequate metho ds. This is the purp ose of ongoing w ork. 7 A kno wledgemen ts It is our pleasure to thank Mi hael Baak e and Natali Zin t for ritially reading the man usript, and Hugo v an den Berg and F rank den Hollander for helpful disussions. This w ork w as supp orted 25 b y DF G-F OR 498 (Dut h-German Bilateral Resear h Group on Mathematis of Random Spatial Mo dels in Ph ysis and Biology) and the NR W In ternational Graduate S ho ol of Bioinformatis and Genome Resear h at Bielefeld Univ ersit y . 26 onstitutiv e(horiz.)-v ariable(v ert.) 7 73 140 207 4 45 87 129 foreign(horiz.)-v ariable(v ert.) 7 76 145 214 4 45 87 129 foreign(horiz.)-onstitutiv e(v ert.) 7 76 145 214 7 87 168 247 13 142 270 399 3 42 82 121 9 100 191 281 3 42 82 121 9 100 191 281 13 167 322 476 16 171 328 485 3 44 85 125 9 101 193 285 3 49 85 125 9 101 193 285 16 203 391 578 18 208 385 569 3 44 85 125 9 101 193 284 4 44 85 125 9 101 193 285 18 238 459 679 31 344 657 970 4 49 88 130 9 101 193 285 4 49 88 130 9 101 193 285 31 407 782 1158 43 478 912 1347 3 44 85 126 9 101 193 285 3 44 85 126 9 101 193 285 43 564 1086 1607 Fig. 13 P airwise join t frequenies of the total onstitutiv e, v ariable, and foreign stim ulation rates, for those samples with G ( z ( f ) ) > g act in the basi mo del (3) (with z ( f ) = 1000 ). Greysales orresp ond to n um b er of samples falling in to 2D-in terv als dened b y total stim ulation rates of pairs of an tigen t yp es. Ro ws (from top to b ottom): g act = 100 , 250 , 350 , 500 , 750 , 1000 ; olumns (from left to righ t): onstitutiv e (horizon tal) v ariable (v ertial); foreign (horizon tal) v ariable (v ertial); foreign (horizon tal) onstitu- tiv e (v ertial). Ligh ter shading orresp onds to higher frequenies. 27 0 0 . 3 0 . 6 0 . 9 300 600 900 g act z ( f ) = 500 z ( f ) = 1000 z ( f ) = 1500 z ( f ) = 2000 Fig. 14 F ration of samples whose self-omp onen t alone is ab o v e threshold, among those that rea h the threshold in the presene of z ( f ) foreign moleules, for v arious z ( f ) (i.e., IS sim ulation of the probabilit y in Eq. (39 )). Sample size is 10000 for ea h g act v alue onsidered. 28 Referenes 1. Arstila, T., Casrouge, A., Baron, V., Ev en, J., Kannelop oulos, J., K ourilsky , P .: A diret estimate of the h uman αβ T ell reeptor div ersit y . Siene 286 , 958961 (1999). 2. Asm ussen, S.: Applied Probabilit y and Queues. 2nd ed., Springer, New Y ork (2003). 3. Billingsley , P .: Probabilit y and Measure. 3rd ed., Wiley , New Y ork (1995). 4. Boro vsky , Z., Mishan-Eisen b erg, G., Y aniv, E., Ra hmilewitz, J.: Serial triggering of T ell reeptors results in inremen tal aum ulation of signaling in termediates. J. Biol. Chem. 277 , 2152921536 (2002). 5. Bu klew, J.A.: In tro dution to Rare Ev en t Sim ulation. Springer, New Y ork (2004). 6. Da vis, S.J., Ik emizu, S., Ev ans, E.J., F ugger, L., Bakk er, T.R., v an der Merw e, P .A.: The nature of moleular reognition b y T ells. Nat. Imm unol. 4 , 217224 (2003). 7. Dem b o, A., Zeitouni, O.: Large Deviations Te hniques and Appliations. Springer, New Y ork (1998). 8. den Hollander, F.: Large Deviations. AMS, Pro videne, RI (2000). 9. Diek er, A., Mandjes, M.: On asymptotially eien t sim ulation of large deviation probabilities. A dv. Appl. Prob. 37 , 539552 (2005). 10. Dushek, O., Co om bs, D.: Analysis of serial engagemen t and p eptide-MHC transp ort in T ell reeptor miro lusters. Bioph ys. J. 94 , 34473460 (2008). 11. Georgii, H.O.: Sto hastis. de Gruyter, Berlin (2008). 12. Gonzalez, P .A., Carreno, L.J., Co om bs, D., Mora, J.E., P almieri, E., Goldstein, B., Nathenson, S.G., Kalergis, A.M.: T-ell reeptor binding kinetis required for T ell ativ ation dep end on the densit y of ognate ligand on the an tigen-presen ting ell. Pro . Natl. A ad. Si. U.S.A 102 , 48244829 (2005). 13. Hla v aek, W.S., Redondo, A., W ofsy , C., Goldstein, B.: Kineti pro ofreading in reeptor-mediated transdution of ellular signals: reeptor aggregation, partially ativ ated reeptors, and ytosoli mes- sengers. Bull. Math. Biol. 64 , 887911 (2002). 14. Hun t, D. F., Henderson, R.A., Shabano witz, J., Sak agu hi, K., Mi hel, H., Sevilir, N., Co x, A.L., App ella, E., Engelhard, V.H.: Charaterization of p eptides b ound to the lass I MHC moleule HLA- A2.1 b y mass sp etrometry . Siene 255 (1992), 12611263. 15. Kalergis, A.M., Bou heron, N., Douey , M.A., P almieri, E., Go y arts, E.C., V egh, Z., Lues her, I.F., Nathenson, S.G.: Eien t T ell ativ ation requires an optimal dw ell-time of in teration b et w een the TCR and the pMHC omplex. Nat. Imm unol. 2 , 229234 (2001). 16. Kronmal, R.A., P eterson, A.J.: On the alias metho d for generating random v ariables from a disrete distribution. Amer. Stat. 33 , 214218 (1979). 17. Lanet, D., Sado vsky , E., Seidelmann, E.: Probabilit y mo del for moleular reognition in biologial reeptor rep ertoires: Signiane to the olfatory system. Pro . Natl. A ad. Si. U.S.A. 90 , 37153719 (1993). 18. Lord, G.M., Le hler, R.I., George, A.J.: A kineti dieren tiation mo del for the ation of altered TCR ligands. Imm unol. T o da y 20 , 3339 (1999). 19. Madras, N.: Letures on Mon te-Carlo Metho ds. AMS, Pro videne, RI (2002). 20. Mason, D.: A v ery high lev el of rossreativit y is an essen tial feature of the T-ell reeptor. Im- m unol. T o da y 19 , 395404 (1998). 21. MKeithan, T.W.: Kineti pro ofreading in T-ell reeptor signal transdution. Pro . Natl. A ad. Si. U.S.A. 92 , 50425046 (1995). 22. Rabino witz, J.D., Beeson, C., W ulng, C., T ate, K., Allen, P .M., Da vis, M.M., MConnell, H.M.: Altered T-ell reeptor ligands trigger a subset of early T ell signals. Imm unit y 5 , 125135 (1996). 23. Rosen w ald, S., Kafri, R., Lanet, D.: T est of a statistial mo del for moleular reognition in biologial rep ertoires. J. Theor. Biol. 216 , 327336 (2002). 24. Ross, S.M.: Sim ulation. A ademi Press (2002). 25. Rothen b erg, E.V.: Ho w T-ells oun t. Siene 273 , 7880 (1996). 26. Sado wsky , J.S., Bu klew, J.A.: On large deviations theory and asymptotially eien t Mon te Carlo estimation. IEEE TIT 36 , 579588 (1990). 27. Sousa, J., Carneiro, J.: A mathematial analysis of TCR serial triggering and do wn-regulation. Eur. J. Imm unol. 30 , 32193227 (2000). 28. Stev ano ví, S., S hild, H.: Quan titativ e asp ets of T ell ativ ation p eptide generation and editing b y MHC lass I moleule. Seminars Imm unol. 11 (1999), 375384. 29. Utzn y , C., Co om bs, D., Muller, S., V alitutti, S.: Analysis of p eptide/MHC-indued TCR do wnregula- tion: deiphering the triggering kinetis. Cell Bio hem. Bioph ys. 46 , 101111 (2006). 30. V alitutti, S., Lanza v e hia, A.: Serial triggering of TCRs: a basis for the sensitivit y and sp eiit y of an tigen reognition. Imm unol. T o da y 18 , 299304 (1997). 31. V alitutti, S., Muller, S., Cella, M., P ado v an, E., Lanza v e hia, A.: Serial triggering of man y T-ell reeptors b y a few p eptide-MHC omplexes. Nature 375 , 148151 (1995). 32. v an den Berg, H.A., Molina-P arís, C.: Th ymi presen tation of autoan tigens and the eieny of negativ e seletion. J. Theor. Med. 5 , 122 (2003). 33. v an den Berg, H.A., Rand, D.A., Burroughs, N.J.: A reliable and safe T-ell rep ertoire based on lo w-anit y T-ell reeptors. J. Theor. Biol. 209 , 465486 (2001). 34. v an den Berg, H.A., Rand, D.A.: An tigen presen tation on MHC moleules as a div ersit y lter that enhanes imm une eay . J. Theor. Biol. 224 , 249267 (2003). 29 35. v an den Berg, H.A., Rand, D.A.: Quan titativ e theory of T-ell resp onsiv eness. Imm unol. Rev. 216 , 8192 (2007). 36. Viola, A., Lanza v e hia, A.: T-ell ativ ation determined b y T-ell reeptor n um b er and tunable thresholds. Siene 273 , 104106 (1996). 37. Zin t, N., Baak e, E., den Hollander, F.: Ho w T-ells use large deviations to reognize foreign an tigens. J. Math. Biol. 57 , 841861 (2008).
Original Paper
Loading high-quality paper...
Comments & Academic Discussion
Loading comments...
Leave a Comment