## Abstract

Asthma is a complex chronic inflammatory disease of the small airways that has dramatically increased in prevalence in industrialized countries during the last decades. Risk factors for adult asthma have been related to the complex array of gene-environment interactions and exposure of the immune system to allergens in early childhood. In genetically predisposed subjects, continuous exposure to environmental agents such as allergens or infections can lead to recurrent airway symptoms characterized by recurrent episodes of airway inflammation and bronchoconstriction with clinical symptoms of cough, dyspnea, or wheezing. In this study, we report that the longterm temporal dynamics of recurrent airway symptoms in a population of unselected infants display a complex intermittent pattern and that the distribution of interepisode intervals follows a power law. We interpret the data by using a model of the dynamics of attack episodes in which an attack is triggered by an avalanche of airway constrictions. We map the dynamics of this model to the known problem of a random walk in the presence of an absorbing boundary in which the walker corresponds to the fluctuations in contractile state of airway smooth muscle cells. These findings may provide new insight into the mechanisms of otherwise unexplained symptom episodes.

- asthma
- instabilities
- distribution
- avalanche
- power law

asthma is a complex chronic inflammatory disease of the small airways that has dramatically increased in prevalence in the industrialized countries during the last decades (4, 31). Risk factors for adult asthma have been related to the complex array of gene-environment interactions and exposure of the immune system to allergens in early childhood (22, 26). In genetically predisposed subjects, continuous exposure to environmental agents such as allergens or infections can lead to recurrent airway symptoms, which are characterized by recurrent episodes of airway inflammation and bronchoconstriction with clinical symptoms of cough, dyspnea, or wheezing (5, 10). The dynamics of the interactions between the immune system and the environment during early infancy can play a critical role in the subsequent development of asthma (26). Epidemiological studies have led to the hygiene hypothesis that Western lifestyle, improved public health and the use of antibiotics in the industrialized countries have reduced the incidence of infections and other environmental exposures of the immune system that would otherwise protect against asthma (16, 18, 27). However, little is known of how the sequence of environmental exposures during early life affects the later development of the immune system (6).

Airway diseases such as asthma and other wheezing disorders with recurrent symptoms may be similar to “dynamic diseases” (2), which can be described by nonlinear dynamics. Two essential facets of nonlinear dynamics are memory and nonlinearity (12). Memory of the immune and inflammatory systems may be apparent in persistent bronchial hyperreactivity after a single viral (9) or allergen challenge of asthmatic patients (3) or persistent symptoms after stimulus removal in occupational asthma (14). Nonlinearity often manifests itself as a threshold mechanism in response to environmental stimuli. When the input to such a dynamic system fluctuates randomly, the response of the system can display complex intermittent behavior (23). If the only observable variable is the output of the system, then predictions of future outputs may not be reliable and can only be given statistically. Indeed, there is often no apparent association between the severity of clinical symptoms and the timing and strength of environmental stimuli, as for example in fatal asthma attacks characterized by a complete collapse of the system (20). Rational therapy would thus certainly benefit from understanding the temporal pattern of gene-environment interactions, but to date predictions of disease incidence or exacerbations have had limited success (20, 34).

The aim of the present study was to characterize the statistical properties of the temporal dynamics of recurrent airway symptoms in early life and to interpret these statistical properties in terms of a model that incorporates the interactions among airway mechanics, cellular function, and the fluctuations in environmental stimuli.

## METHODS

Characterization of recurrent airway symptoms during the first year of life associated with wheezing disorders and asthma has been difficult since it requires tracking the sequence of exacerbation episodes for an extended period of time on an individual basis. On the basis of clinical grounds, wheezing disorders in infants with mainly episodic viral wheeze, cough, or dyspnea, which are often transient in the first years of life, are difficult to distinguish from chronic asthma with recurrent lower airway symptoms and persistence into school age. Clinical symptoms during an episode are similar in both phenotypes since viral infections are the most frequent triggers causing airway inflammation, mucus secretion, or reversible bronchoconstriction. In an unselected birth cohort, both phenotypes will be present. A further difficulty arises from the assessment of clinical symptoms by parental judgement. In general, parents weigh symptoms differently, which can introduce subjective errors in the collected data. To minimize these difficulties, established respiratory symptom scores (25) were assessed by the same trained research assistants by using standardized weekly interviews in a prospective birth cohort study of 78 unselected Caucasian infants. An episode of lower respiratory disease was defined and scored (0 if no symptom, 1 or higher in the presence of symptoms) based on symptoms of cough, wheeze, or dyspnea if they lasted for at least 2 consecutive days. All infants entered the study at birth with no birth-related respiratory problems and stayed in the study for 1 yr involving all four seasons. Although the infants likely had different genetic backgrounds, the subject population was relatively homogeneous in terms of study starting conditions and environmental exposures.

The study was approved by the Institutional Review Boards of both the University of Berne and Boston University. In each case, written consent was obtained from the accompanying parent.

## RESULTS

*Experimental data.* We analyzed the periods between two consecutive episodes, interepisode intervals (IEI), defined as the number of days between two nonzero scores (Fig. 1*A*). On average, infants had seven (ranging from 1 to 23) episodes in a year, providing altogether 566 IEI values. This is illustrated by the representative examples of IEI time series in two infants in Fig. 1, *B* and *C*. The infant in Fig. 1*B* had 7 IEI values representing the average number of IEI per year, whereas the infant in Fig. 1*C* had the longest sequence among the infants with 23 IEI values in the first year of life. Notice the significant fluctuations and intermittent behavior in both cases with many small and a few large IEI.

On the basis of the homogeneity of the cohort, we pooled the 566 IEI values from all subjects (Fig. 2*A*) to construct a single distribution of IEI representative of this infant population (Fig. 2*B*). This distribution shows a peak around IEI = 10 days. It is unclear whether, for IEI of <10 days, two episodes belong to a single episode or whether they are two separate episodes. Thus we only consider the data for IEI of >10 days for which the distribution follows a power law (straight line on log-log graph) over nearly 1.5 decades with an exponent (negative slope) α = 1.51 ± 0.14.

We examined the possibility that some of the scores in Fig. 1*A* provided by the parents were incorrect. First, we randomly selected 10% of the IEI intervals. Each interval was then split randomly into two subintervals to mimic the case when a score of 1 was reported by the parents, although there were no real clinical symptoms. The effect was to increase the likelihood of small IEI, but the tail and α of the distribution remained the same. Next, we considered the case when a score of 0 was reported, although there were noticeable clinical symptoms. In this case, we randomly selected 10% of the IEI intervals and coalesced two neighboring intervals, which had a small effect on the tail but did not change α appreciably. Finally, we combined the two scenarios and found that even if 20% of the intervals were not correctly identified, the error in exponent would be <5%. We conclude that the distribution of IEI is robust against small subjective errors during the interviews.

*Model development.* To interpret the power law distribution of IEI, we developed a computer model of how continuous exposure to environmental agents can lead to intermittent fluctuations in clinical symptoms (Fig. 3). The lung is modeled as a three-dimensional bifurcating structure similar to the tracheobronchial tree (13). To account for the dynamic interplay between mediators, remodeling, and inflammation (11), we incorporate into each segment a functional unit related to the inflammatory and remodeling state of the airway and the contractile properties of the smooth muscle in the airway wall. When environmental agents are inhaled, they are distributed randomly in the airways. In the model, we assume that each segment of the airway tree receives a fraction *s* of the total amount of the inhaled agents (*E*), and these agents can reduce the cross-sectional area of the airway by contracting the smooth muscle and/or secretion of mucus, which occurs in the model when *s* > 0. The effects of decreasing external stimuli, smooth muscle relaxation, and elimination of mucus, cellular debris, and inflammation due to internal cellular repair (30) are modeled by *s* < 0. Thus the time series of *s* is a zero mean uniformly distributed white noise.

The effects of various agents can continuously build up at the cellular level and involve remodeling, inflammation, and changes in the contractile state of smooth muscle (7, 8, 17, 19). These processes occur on rather different time scales. Remodeling is a slow adaptive process that may take months and years. Inflammation is much faster and can develop with minutes or hours, whereas changes in contractile state are rapid, nearly instantaneous. In the model, we take these effects into account by assigning memory to the functional units in the sense that the instantaneous state of each unit depends not only on the present value of *s* but on the time history of the stimulus. We quantitatively model this phenomenon by storing the cumulative sum *S* = Σ*s*(*k*) of the past stimuli for each segment where the sequence *s*(*1*), *s*(*2*),... *s*(*M*) constitutes the time series of stimuli, with *s*(*M*) being the most recent value of *s*. After a build up of the effects of the stimuli, a small additional stimulus is sufficient to generate a disproportionate cellular response, which we model by setting up a threshold *T*. The threshold *T* in the model is related to the mechanical properties of the airway wall, the fragility of the bronchial epithelium, and the presence of mucus, inflammation, etc. If the next stimulus *s*(*M* + 1) increases *S* for a given segment such that *S* > *T*, the smooth muscle of that segment is maximally contracted. An airway with maximally contracted smooth muscle behaves as a bistable system, i.e., it is either fully open or nearly collapsed (1). Whether an airway jumps to the open or closed state after maximal contraction depends on several factors, including the elastic support (or tethering) of the surrounding parenchyma. Because our model is three dimensional, we can also take into account the spatial interactions among adjacent segments. When an airway (labeled *i*) collapses (*S*_{i} > *T*), it pulls a neighboring airway *j* open, which in turn releases parenchymal tethering around itself and potentially allows another spatially adjacent airway (labeled *k)* to collapse. Whether the collapse of airway *i* triggers the collapse of airway *k* depends on the state *S*_{k} of airway *k*. We model the effects of reduced tethering on the state of airway *k* by adding an additional random amount of *s* to *S*_{k}. If the new value of *S*_{k} satisfies the condition *S*_{k} > *T*, then airway *k* also collapses, which in turn may also trigger the opening of another airway. The process of sequential collapse-open-collapse in neighboring airways can propagate spatially, and hence maximal muscle contraction in an airway can generate a cascade or avalanche of bistabilities and constrictions. If such an avalanche propagates through a large portion of the tree, a substantial decrease in airway lumen and an increase in airway resistance occur. To mimic the scores in the measured data, we assign a score of 1 to a particular iteration if an avalanche was triggered during that step. After an avalanche, we allow the system to relax by erasing the memory of each unit (*S*_{i} = 0 for all *i*). Next, loading the system with environmental agents starts again. A schematic algorithm of the model is summarized in Fig. 3.

*Numerical simulations.* We simulate these processes numerically in a tree model and find that the model displays an intricate sequence of avalanches. We define IEI as the interval, or number of iterations, between two consecutive avalanches. Because the model incorporates many essential features of nonlinear dynamic systems including memory, threshold phenomena, spatial interactions, and avalanches, the simulated times series of IEI show complex intermittent behavior (Fig. 4*A*) similar to the data in Figs. 1 and 2. Additionally, the distribution of IEI from the model is a single power law with α = 1.54 ± 0.05, in agreement with the experimental data (Fig. 4*B*).

## DISCUSSION

*Relation to random walk.* In general, systems with avalanche dynamics often display intermittent behavior with the distribution of waiting times between avalanches following a power law (24, 28, 29, 33). The power law distribution in our model arises from the following ingredients. The continuous competition between the buildup (increasing *S* mimicking remodeling, inflammation, and muscle contraction) and the relaxation (decreasing *S* mimicking cellular repair) leads to fluctuations in the state of the functional unit that slowly approaches the threshold *T* (maximal contraction). The fluctuations in muscle state are then terminated by an avalanche-like spatial propagation of bistabilities causing the intermittent behavior of IEI. These fluctuations can be considered as a random walk. A small amount of muscle contraction or relaxation corresponds to a walker taking a unit step in the positive or negative direction, respectively. The contractile state of the muscle (*S* in the numerical model) is analogous to the position or end point of the walk, which is the cumulative sum of the unit steps. The smooth muscle of each airway segment can be considered as a separate random walker. When the smooth muscle reaches maximum contraction, an avalanche is triggered and the system is reset, which is equivalent to resetting the walker to the origin when it reaches a threshold, i.e., the cumulative sum of the unit steps reaches a given positive number. The resetting itself is similar to the problem of a random walk with an absorbing boundary. In this case, the walker steps up and down with equal probability at every iteration, but once it reaches a preset distance from the origin, it cannot return, hence the name absorbing boundary. The survival time of the walker is the number of steps required to reach the absorbing boundary, which is analogous to the number of iterations to reach an avalanche, or the IEI, in our numerical model. The distribution of IEI is thus equivalent to the distribution *P* of the survival time τ, which is known to have an analytical solution (32) 1 where *S*_{0} is the initial distance of the walker from the absorbing boundary, and *D* is the diffusion coefficient. The distribution of survival times is a power law with an exponent of 1.5 in agreement with our data and numerical simulations. Thus the distribution of IEI during a 1-yr follow-up of recurrent airway symptoms is consistent with statistical properties of a cellular random walk.

*Physiological implications of the modeling.* The possible implications of *Eq. 1* are intriguing. Asthmatic subjects have clinical symptoms substantially more often than healthy subjects, and hence one would expect that the probability of a given IEI is smaller for asthmatic than for normal subjects. Fluctuations in airway resistance have been shown to be due to the fluctuations in smooth muscle contractile state, and the fluctuations in muscle length have been argued to be faster and exaggerated in asthmatic compared with healthy individuals (21). Thus, to the extent that *Eq. 1* holds for asthmatic subjects, one may speculate that the decreased likelihood of long IEI in asthmatic subjects may be due to a smaller initial distance *S*_{0} and/or an increased diffusivity *D*. Changes in these parameters may represent remodeling of the airway wall or hypertrophy of the muscle (17), either of which facilitates the collapse of an airway, or the reactivity of the airways due to enhanced inflammatory response, increased smooth muscle shortening velocity (21), cytoskeletal remodeling (8), or perhaps a stiff and very contractile smooth muscle frozen in a so-called latch state (7). The clinical consequences may be expressed as bronchial hyperreactivity in asthmatic subjects. It is of course also possible that *Eq. 1* simply does not hold for asthmatic subjects or that the exponent in *Eq. 1* increases in asthmatic subjects. In either case, such a mechanistic approach to the statistical fluctuations in the observed IEI distribution should allow for a better prediction of future attack episodes on an individual basis than previous statistical methods that are based on epidemiological studies.

It is important to emphasize that the good correspondence between the measured and simulated IEI distributions in Figs. 2*B* and 4*B* does not necessarily imply that the experimentally obtained distribution is determined by the same mechanisms that are incorporated into our computer model. Additionally, we also note that the quantitative agreement between our data and the random walk model does not directly prove that the temporal dynamics of recurrent airway symptoms are a consequence of a cellular random walk. Nevertheless, the model presented here represents a first attempt to link the interactions among cellular function, environmental stimuli, and airway tree structure to the statistical distribution of clinically measurable symptom episodes in a population of subjects. We believe that such modeling efforts will help bridge the gap between basic cell/molecular biology and the clinical sciences, which, in recent years, appear to have been undergoing divorce (15).

In conclusion, in this study, we have introduced a new statistical approach to analyze the temporal dynamics of recurrent airway symptoms. The statistical analysis of IEI may be a useful technique to quantify fluctuations and exacerbations of various other recurrent diseases. Our findings may also provide new insight into the mechanisms of otherwise unexplained symptom episodes and have implications for rational therapy or preventative measures of recurrent inflammatory diseases and the design of outcome measures in future epidemiological studies.

## DISCLOSURES

This study was supported by the Swiss National Science Foundation and National Science Foundation Grant BES 0114538.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2003 the American Physiological Society