Chaotic dynamics and diffusion in a piecewise linear equation

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 CHAOS 25, 033103 (2015) Chaotic dynamics and diffusion in a piecewise linear equation Pabel Shahrear,1,a) Leon Glass,2,b) and Rod Edwards3,c) 1 Department of Mathematics, Shah Jalal University of Science and Technology, Sylhet–3114, Bangladesh Department of Physiology, 3655 Promenade Sir William Osler, McGill University, Montreal, Quebec H3G 1Y6, Canada 3 Department of Mathematics and Statistics, University of Victoria, P.O. Box 1700 STN CSC, Victoria, British Columbia V8W 2Y2, Canada 2 (Received 7 December 2014; accepted 11 February 2015; published online 10 March 2015) Genetic interactions are often modeled by logical networks in which time is discrete and all gene activity states update simultaneously. However, there is no synchronizing clock in organisms. An alternative model assumes that the logical network is preserved and plays a key role in driving the dynamics in piecewise nonlinear differential equations. We examine dynamics in a particular 4dimensional equation of this class. In the equation, two of the variables form a negative feedback loop that drives a second negative feedback loop. By modifying the original equations by eliminating exponential decay, we generate a modified system that is amenable to detailed analysis. In the modified system, we can determine in detail the Poincar!e (return) map on a cross section to the flow. By analyzing the eigenvalues of the map for the different trajectories, we are able to show that except for a set of measure 0, the flow must necessarily have an eigenvalue greater than 1 and hence there is sensitive dependence on initial conditions. Further, there is an irregular oscillation whose amplitude is described by a diffusive process that is well-modeled by the Irwin-Hall distribution. There is a large class of other piecewise-linear networks that might be analyzed using similar methods. The analysis gives insight into possible origins of chaotic dynamics in periodically C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4913417] forced dynamical systems. V Genetic control networks are fundamental to the proper functioning of organisms. Simplified mathematical models of these networks assume that genes are turned on and off abruptly as molecules that control gene expression cross critical thresholds. The resulting class of differential equations can show a variety of different dynamic behaviors including steady states, oscillations, and chaos. We have been interested in developing methods to analyze the resulting class of differential equations. The current paper treats one particular example that has unusual dynamics in which trajectories assume irregular spike-like and sawtooth trajectories. By simplifying the original equations, we are able to understand the source of the irregular trajectories. In particular, two of the variables have oscillatory dynamics with a variable amplitude. We show that the amplitude of the oscillation fluctuates by a diffusive process. These results give further insight into the source of irregular dynamics in complex nonlinear networks such as those that control key biological functions. I. INTRODUCTION Gene expression in organisms is controlled by specific protein molecules called transcription factors that regulate their activity. Depending on the gene, the transcription factors present, and the control function, genes may either be a) [email protected] [email protected] c) [email protected] b) 1054-1500/2015/25(3)/033103/15/$30.00 “switched on” or “switched off”. Over 40 years ago, Kauffman proposed that these gene networks could be modeled by Boolean networks with discrete time and synchronous updating.1 However, since organisms do not possess clocks that would enable synchronous updating, there has been interest in the study of a class of piecewise-linear differential equations in which the logical structure of the Boolean network is preserved.2 In particular, such equations can display nonlinear behavior including limit cycle oscillation3–6 and chaotic dynamics.7–9 In unpublished work, one of us (R.E.) carried out a search of randomly constructed 4 dimensional equations with no self-input. He found a small number of networks that had bizarre dynamics manifested by spike-like and square wave trajectories when the phase space portrait was projected down to two dimensions. Since such dynamics had not been described previously, we have initiated studies of the equations. In a recent paper, we studied a particular equation as a parameter that controls the steepness of the control functions was varied.10 The bizarre geometries were preserved for continuous versions of the piecewise linear equations provided the control functions were sufficiently steep. But, this work did not provide insight into the dynamical behavior. In the following, we study another such 4-dimensional system with bizarre dynamics. By modifying the original equations, we obtain a simpler system that enables a detailed understanding of the dynamics. In the modified system, two of the variables form a negative feedback loop and act as an input to the other two variables. The other two variables display a complex oscillation with fluctuating amplitude. A 25, 033103-1 C 2015 AIP Publishing LLC V This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-2 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) single cycle of the oscillation can be represented by 16 different symbolic sequences. We can determine regions of state space associated with each sequence. We also determine the eigenvalues associated with each sequence and show that there must be one eigenvalue greater than 1 so that the map is expanding. We have strong analytical and numerical evidence for chaotic dynamics in the flow. The chaotic dynamics are manifest by a fluctuating amplitude of the driven oscillator. The statistical properties of this amplitude fluctuation are well described analytically by the Irwin-Hall distribution, which characterizes the distribution of random sums of numbers from the unit interval.11,12 In Sec. II, we present the original 4-dimensional equations and some numerical results concerning the dynamics. These results led us to reformulate the equations. In Sec. III, we present the reformulated system and show how to compute the dynamics analytically. In Sec. IV, we compute how the paths through phase space depend on the intitial conditions. Based on the analyses in Secs. III and IV, we find strong evidence that the dynamics are chaotic. We also show that the fluctuations of amplitude of the oscillations are well described as a diffusion process using the Irwin-Hall distribution.11,12 In Sec. V, we partition the space of initial conditions (Poincar!e section) into regions from which different paths are followed and on which different return maps apply. In the Appendix, we provide additional information concerning the trajectories, a proof that the model in Secs. III and IV is senstive to initial conditions, and analysis of small amplitude trajetories. II. A FOUR-DIMENSIONAL PIECEWISE LINEAR EQUATION The following set of equations was found in studies of randomly generated 4-dimensional piecewise linear models of gene networks like those initially studied by Glass and coworkers.2,3 dx1 dt dx2 dt dx3 dt dx4 dt TABLE I. Truth Table for Eq. (1). X1 X2 X3 X4(t) 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 X1 X2 X3 X4(t þ 1) 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 1 0 0 0 1 0 0 1 1 0 1 1 1 0 1 1 0 1 0 1 0 1 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 Boolean state to the subsequent state in the truth table. As shown previously,3 provided there is no self-input, each edge in the resulting state transition diagram has a unique orientation. Figure 1 shows the associated state transition diagram on the 4-cube for this network. The state transition diagram reveals a complex structure with no “sinks” (vertices in which all edges are directed inwards), and no attracting cycles. Both these structures are important since their presence allows one to infer stable steady states or stable limit cycles, respectively, in the associated differential equations.3 However, the presence of multiple return cycles to the vertices shows that the system satisfies a necessary condition for chaotic dynamics.9 ¼ 2X 4 " 1 " x1 ; ! " ¼ 2 X 3 X1 þ X1 X3 X4 þ X 1 X 3 X4 " 1 " x2 ; ! " ¼ 2 X 1 X 4 þ X1 X2 " 1 " x3 ; (1) ¼ 2X1 " 1 " x4 ; where xi are continuous variables, Xi are logical variables defined by Xi ¼ H(xi), where H(x) is the Heaviside step function, X"i ¼ 1 " Xi for each i ¼ 1,…,4. As described previously,3 these equations are based on an underlying logical structure that can be expressed as a truth table. The truth table for these equations is presented in Table I. One representation of the network is provided by mapping the logical structure of the network onto a Boolean 4cube. Each vertex of the hypercube represents an orthant of phase space. Directed edges on the hypercube show allowed flows in the four-dimensional phase space. The outdegree of each vertex is equal to the Hamming distance from a FIG. 1. A hypercube representation of the network based on the truth table. Each vertex represents a logical state, and the edges between the vertices represent the allowed transitions in the associated differential equation. This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-3 Shahrear, Glass, and Edwards Numerical integration of these equations shows complex dynamics. Figure 2 shows a projection of the dynamics in each of the six two-dimensional subspaces. Figure 3 shows segments of the trajectories over two different time intervals. These results provide a clear conceptual picture of the structure of the network. The amplitude of the oscillation of x1, x4 decays gradually as seen in Figure 3. The projection in the x2, x3 subspace shows a complex oscillatory trajectory. Understanding this x2, x3 projection is the main focus of the following. In order to better appreciate this structure, it is useful to rewrite the truth tables. Variables 1 and 4 form a simple negative feedback loop (see Table II). As described and proved in earlier work,13 variables 1 and 4 cycle and asymptotically approach the origin. Because of the slow rate of decay, the decay is not obvious in the x1, x4 projection in Fig. 2. As a consequence of the decay of the amplitude, the period of the oscillation continually decreases, and this acts as an input to x2 and x3. The dynamics of variables 2 and 3 are more complex. To help understand these dynamics we construct the logical truth tables that control them in Tables III and IV. Based on the truth tables for x2 and x3 there is a negative feedback loop in which x2 activates x3, and x3 inhibits x2. However, there is also a forcing from the cyclic x1, x4 oscillation. The result of this forcing is that when x2 is high, the trajectory of x3 consists of 3 steps up and one step down; and when x2 is low, the trajectory of x3 consists of 3 steps down and one step up. Conversely, when x3 is high, the trajectory of x2 consists of 3 steps down and one step up; and when x3 is low, the trajectory of x2 consists of 3 steps up and one step down. The phase space in the x2, x3 projection shows a “square” geometry in which there are spikes or sawteeth along some of the edges. To avoid having these short term fluctuations confound the analysis we define the amplitude of the oscillation as the minimum value of jx2j þ jx3j for Chaos 25, 033103 (2015) three consecutive steps. One measure of the long time dynamics is the amplitude of the x2, x3 oscillation. A numerical appreciation of the fluctuation of amplitude is shown in Figure 4. Although there is a gradual decline in the amplitude, the decrease is not monotonic. III. A MODIFIED EQUATION WITH PERIODIC FORCING To analyze the dynamics, we consider a simplified model in which x1, x4 oscillate with fixed period and amplitude. For x1, x4, we assume that the period of the oscillation is 4 and that each element oscillates in the range 61. Further, for x2 and x3, we eliminate the exponential decay terms so that each variable is either increasing with a slope of þ1 or decreasing with slope of "1. Thus, the x1, x4 oscillation is a periodic input that forces the x2, x3 subsystem. The modified equations derived from Eq. (1) are dx1 dt dx2 dt dx3 dt dx4 dt ¼ 2X 4 " 1; ! " ¼ 2 X 3 X1 þ X1 X3 X4 þ X 1 X 3 X4 " 1; ! " ¼ 2 X 1 X 4 þ X1 X2 " 1; (2) ¼ 2X1 " 1; where we assume an initial condition of (x1,x2,x3,x4) with the criterion jx1j þ jx4j ¼ 1. Our goal is to achieve a detailed understanding of the dynamics of this system. We discuss the relationship of this reduced system to the original one in Sec. VI. In most of this paper, we assume a “large amplitude limit” in which the amplitude of the oscillation is greater than 2. By analyzing algebraically the flows in the x2, x3-plane, we are able to show that for most initial conditions, there is sensitive dependence on the initial conditions so that two initial conditions that are close together diverge over time. The FIG. 2. Projections of Eq. (1) in the 6 subspaces defined by taking each pair of two variables. The projections are plotted for 300 threshold crossing following a transient of 9000 threshold crossings from an initial condition of: x1 ¼ 0.500, x2 ¼ 0.036, x3 ¼ "0.057, x4 ¼ 0.500. This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-4 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) FIG. 3. Representative traces of the dynamics in Eq. (1). x1 (dark blue), x2 (green), x3 (red), x4 (light blue). TABLE II. Truth table for a negative feedback between x1 and x4. X1 X4(t) 0 0 1 1 0 1 0 1 X1 X4(t þ 1) 1 0 1 0 0 0 1 1 basic picture that emerges is that the amplitude of the oscillation in the x2, x3-plane undergoes a “diffusive” process. In each orbit in the x2, x3-plane, there is one increment in amplitude in the interval [0, 2] and a second increment in the interval ["2, 0], and these increments have the appearance of being selected randomly from a uniform distribution. We show the two-dimensional projections of a trajectory in our modified system onto the x1, x4-plane and the x2, x3plane in Figure 5. To study the dynamics, we consider the trajectory in the x2, x3-plane. It has a basic “square” TABLE III. Logical structure of X2 which depends on X3X4X1. X3 X4 X1(t) 0 0 0 0 1 1 1 1 0 0 1 1 0 0 1 1 0 1 0 1 0 1 0 1 X2(t þ 1) 0 1 1 1 0 0 0 1 TABLE IV. Logical structure of X3 which depends on X2X4X1. X2 X4 X1(t) 0 0 0 0 1 1 1 1 0 0 1 1 0 0 1 1 0 1 0 1 0 1 0 1 X3(t þ 1) 1 0 0 0 1 1 0 1 geometry. In the following, we will refer to the “edges” and “corners” of the “square” dropping the quotation marks. The south corner is labeled I, the east corner II, the north corner III, and the west corner IV. Along each edge of the square trajectory, x2 and x3 take 3 steps of unit length in one direction (either up or down) and one step in the other direction. Each step begins and ends at transition points of x1 or x4. The geometry of the trajectory along each edge (either spikes, smooth, or square wave) results from the relative phases of the dynamics of x2 and x3 as determined from the truth tables. In order to fix ideas, in Figure 6, we plot jx2j þ jx3j for a representative part of the trajectory, where we color code the trajectory by the quadrant in the x2, x3-plane. From this plot, we see that the amplitude only changes at corner III (decrease in amplitude from green to black dots) and corner IV (increase in amplitude from black to red dots). Consider the trajectory of x3 along the edge from corner I to corner II. Along this edge, we define a symbolic sequence…, u1,u2,u3,d,u1,u2,u3,d, … where the ui represents steps up and d represents a step down. Crossing of the threshold where x3 ¼ 0 can only first occur during u2 or u3. If the crossing occurs during u2, the next step will also be up so that the threshold will be crossed once. However, if the crossing occurs on step u3, the threshold will be crossed 3 times on u3, d, and u1. Consequently, there are two ways each corner can be traversed: a simple way in which the threshold is crossed once (path A) and a more complicated way in which the threshold is crossed three times (path B). We track the dynamics in the x2, x3-plane. For x2 and x3, we adopt superscripts of the form J, j where J 2 fI; II; III; IVg designates the corner and j 2 fi; og designates whether the value of the variable is at the input or output of a given corner, where these points are defined to be at x1 ¼ 0 represents the value of x3 that is the or x4 ¼ 0. Thus, xIII;i 3 input to corner III. When appropriate, we use vector notation to indicate the input and output coordinates, e.g., xI;i is a colI;i umn vector with two components xI;i 2 and x3 . To express transitions, we also use the ceiling function due (smallest integer $ u) and the floor function buc (largest integer % u). Determination of the transitions at each corner is an arithmetic problem, based on the truth tables. In order to illustrate the dynamics, we consider a transition that occurs This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-5 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) FIG. 4. Fluctuations of the amplitude in Eq. (1). FIG. 5. Two-dimensional projection of a trajectory in the modified equation onto the x1, x4-plane (left panel), and the x2, x3-plane (right panel). at Corner I, starting at an initial value of ð1; "a; "b; 0Þ, where 0 < a < 1 and b > 1 at time t0. From Tables III and IV, we see that in the quadrant where x1 > 0 and x4 > 0, x2 is increasing and x3 is decreasing. Hence, the next transition when variable x2 crosses zero occurs at time t0 þ a when the values of the variables are ð1 " a; 0; "b " a; aÞ. Finally, the next time, a variable equals 0 occurs at time t0 þ 1, when the values of the variables are ð0; 1 " a; 1 " 2a " b; 1Þ. The change of amplitude at the corner is ðð1 " aÞ þð"1 þ 2a þ bÞ "ða þ bÞ ¼ 0. This completes the transition through corner I. On the edge between corners I and II, both x2 and x3 show net increases. Since—at the exit of corner I—x2 >0 and x3 < 0, the time to reach corner II is set by x3. Specifically, we find that I;o I;o xII;i 2 ¼ x2 " dx3 e ; I;o I;o xII;i 3 ¼ x3 " dx3 e: In the Appendix, we summarize the transitions through each of the four corners. Based on this analysis, we compute that the amplitude only changes at corners III and IV. At corner III, calling a ¼ xIII;i 2 , the change in amplitude is "2a on path A, and 2"2a on path B. Both of these give a value between 0 and "2, depending on the value of a. At corner IV, calling b ¼ xIV;i 3 , the change in amplitude is 2b on path A, and 2–2b on path B. Both of these give a value between 0 and 2. IV. SUMMARIES OF TRANSITIONS AND MAPS FIG. 6. jx2j þ jx3j, for a portion of a trajectory in the modified equation. Amplitudes are plotted in colours representing the quadrant in the x2, x3plane: x2 < 0, x3 < 0, red dots; x2 > 0, x3 < 0, blue dots; x2 > 0, x3 > 0, green dots; x2 < 0, x3 > 0, black dots. In the preceding and in the Appendix, we found that based on the values of the variables at the inputs of each corner, we could determine the path at the corner, the input- This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-6 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) TABLE V. Summary of transition paths at each corner. The value in each table is the path taken at the corner, given the path taken at the previous corner and the output point of the previous corner. Corner II path bxI;o 3 c mod 2 Corner I path A B 0 1 B A A B Corner III path bx2II;o c mod 2 Corner II path A B 0 1 A B B A Corner IV path bx3III;o c mod 2 Corner III path A B 0 1 B A A B Corner I path bxIV;o 2 c mod 2 Corner IV path A B 0 1 A B B A output maps for the transitions at each corner, and the maps from the output of one corner to the input of the next corner. By combining these results, we are able to provide a comprehensive algebraic summary of the dynamics in this system. Table V provides a summary of the different paths taken at each corner, given the path taken at the previous corner and the output point of the previous corner. We can also compute the mappings for each of the different circuits from an input of corner I to the next input of I;i corner I, Table VI. Let xI;i 2 ¼ "1 þ r and x3 ¼ "k þ s, where 0 < r < 1, 0 < s < 1 and k > 1 an integer. Note that bxc ¼ x " fxg and dxe ¼ x " fxg þ 1, except when x is an integer. Also, f"xg ¼ 1 " fxg. Based on the above, we can determine the affine maps that take values of r, s at the input to Corner I, on one circuit, to the values of r, s at the input to Corner I at the next circuit. At each corner, we have described two different paths, A and B. Hence, one circuit can be described by a sequence of four letters LIV LIIILIILI, where Lj ¼ A if path A is taken at corner j, Lj ¼ B if path B is taken at corner j, and Lj ¼ x if either path A or B is taken at corner j. The matrices in the mappings for the possible cycles are summarized in Table VII. Consequently, with the exception of the map associated with paths AAxx, the maps have an exponentially expanding direction, because there is an eigenvalue with absolute value jk1j > 1, as well as a contracting direction, with eigenvalue jk2j < 1. The image of a point under one of these maps may of course be in a region where one of the other maps applies. Thus, local expansion and contraction under iteration of the full mapping depends on arbitrary products of the matrices of the four maps. Since each of the four matrices has determinant 1, the determinant of any product matrix (and thus the product of its eigenvalues) is 1. This guarantees that unless jk1j ¼ jk2j ¼ 1 (as is the case for the any powers of the AAxx matrix and only these), we have jk1j > 1 and jk2j < 1, and thus an expanding direction. Proofs of these assertions are in the Appendix. The set of initial points from which we follow paths with the symbolic sequence AAxx indefinitely is shown TABLE VI. Mappings around the cycle. Step I, i I, o II, i II, o III, i III, o, A III, o, B IV, i, AA IV, i, BA IV, i, AB IV, i, BB IV, o, AA IV, o, BA IV, o, AB IV, o, BB I, i, AA I, i, BA I, i, AB I, i, BB x2 "1 þ r r k " r " s þ f2r þ sg k þ 1 " r " s " f2r þ sg 1 " f3r þ 2sg "f3r þ 2sg "f3r þ 2sg "k " 1 þ r þ s þ f2r þ sg " 2f3r þ 2sg "k þ 1 þ r þ s þ f2r þ sg " 2f3r þ 2sg "k " 2 þ r þ s þ 2f3r þ 2sg " f10r þ 7sg "k þ r þ s þ 2f3r þ 2sg " f10r þ 7sg "k þ r þ s þ f2r þ sg " 2f3r þ 2sg "k þ 2 þ r þ s " 3f2r þ sg " 2f3r þ 2sg "k " 1 þ r þ s þ 2f3r þ 2sg " f10r þ 7sg "k " 3 þ r þ s þ 2f3r þ 2sg þ 3f10r þ 7sg "f3r þ 2sg "f11r þ 6sg "f3r þ 2sg f37r þ 26sg " 1 x3 "k þ s "1 " k þ 2r þ s f2r þ sg " 1 f2r þ sg k " r " s þ f3r þ 2sg k þ 1 " r " s þ f3r þ 2sg k þ 1 " r " s " 3f3r þ 2sg f2r þ sg f2r þ sg 1 " f10r þ 7sg 1 " f10r þ 7sg f2r þ sg " 1 f2r þ sg " 1 "f10r þ 7sg "f10r þ 7sg "k " 1 þ r þ s þ 2f2r þ sg " f3r þ 2sg "k þ 1 þ r þ s " 2f2r þ sg " 2f3r þ 2sg þ f11r þ 6sg "k " 1 þ r þ s þ 3f3r þ 2sg " 2f10r þ 7sg "2 " k þ r þ s þ 2f3r þ 2sg þ 2f10r þ 7sg " f37r þ 26sg x1 xI;i 1 2 f0; 1g 1 " xI;i 1 I;o ðxI;o 1 þ dx3 eÞ mod 2 1 " xII;i 1 "ððxII;o þ bxII;o 1 2 cÞ mod 2Þ "1 " x1III;i "1 " x1III;i ðx1III;o þ bx3III;o cÞ mod 2 ðx1III;o þ bx3III;o cÞ mod 2 ðx1III;o þ bx3III;o cÞ mod 2 ðx1III;o þ bx3III;o cÞ mod 2 1 " xIV;i 1 1 " xIV;i 1 1 " xIV;i 1 1 " xIV;i 1 ðxIV;o þ dxIV;o 1 2 eÞ mod 2 ðxIV;o þ dxIV;o 1 2 eÞ mod 2 IV;o ðx1 þ dxIV;o 2 eÞ mod 2 ðxIV;o þ dxIV;o 1 2 eÞ mod 2 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-7 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) TABLE VII. Expansion and contraction factors for the four mappings. Map AAxx BAxx ABxx BBxx Matrix $ $ $ $ "3 "2 2 1 Eigenvalues % "3 "2 "10 "7 "11 "6 2 1 37 26 "10 "7 % % % "1; "1 pffiffiffiffiffi "56 24 pffiffiffiffiffi "56 24 pffiffiffiffiffiffiffiffi 156 224 below to be of measure zero, so, in general, trajectories will be sensitive to initial conditions. However, since the dynamics is not bounded, we cannot assert that the behavior is chaotic. We conjecture that for large amplitudes, the dynamics can be described as follows. There is a continuous circulation around the x2, x3-plane. On each circuit, the amplitude decreases by a random amount in the interval [ " 2, 0] and increases by a random amount in the interval [0, 2]. Thus, there appears to be a “diffusive process” driven by a deterministic dynamical system. The Irwin-Hall distribution11,12 is the distribution of n randomly chosen numbers uniformly distributed in the interval [0, 1]. By the central limit theorem, the sum of many random variables with this distribution is asymptotically Gaussian. For large n, the density is approximated by $ % 1 ðv " lÞ2 ; (3) PðvÞ ¼ pffiffiffiffiffiffi exp " 2r2 r 2p where v is the value of the variable, P(v) is the probability density function, l ¼ n/2 is the mean, and r2 ¼ n=12. For the current situation, based on the preceding analysis, we assume that on each circuit, the amplitude changes by a random number uniformly distributed in ["2, 0] at corner III and in [0, 2] at corner IV. Consequently, for one circuit, the density distribution is a sum of random variables whose distribution is given by 8 1 1 > > þ z; if " 2 % z < 0 > > > <2 4 f ðzÞ ¼ 1 1 " z; if 0 % z < 2 > > > 2 4 > > : 0; otherwise: In this case, for a finite number of steps, N, one has the Irwin-Hall distribution,11,12 in which n ¼ 2N and r2 ¼ n=3. However, as we show below in Sec. V, assuming a uniform IV;i 1 density of the values fxIII;i 2 g and fx3 g, on each circuit 12 of the points will return to the same value after two circuits of the phase space. Consequently, the effective value used in the Irwin-Hall distribution is reduced so that n ¼ 2N–N/3 (notice that since each circuit is equivalent to two steps of the IrwinHall process, return to the same point after two circuits reduces the number of mixing steps of the Irwin-Hall process by 4). Figure 7 compares the theoretically computed distribution based on the N ¼ 30 circuits of the phase space with the theoretically predicted Irwin-Hall distribution ðn ¼ 2 ( 30 "30=3 ¼ 50Þ, starting from an initial amplitude of 5p on the edge between corner IV and corner I for 10 000 initial points. There is close agreement between the predicted distribution and the numerical simulation of the modified equation. Since our process has amplitudes bounded below by 0, the diffusion of amplitudes has a minimum at 0. For this reason, we chose a large amplitude for which this has small influence. The theory provides a good description of the distribution of the fluctuations of amplitude. Small discrepancies may arise from the reflecting boundary where the amplitude is 0, or variation from a uniform invariant density. The preceding analysis and numerical simulations show that assuming an initial large amplitude for the oscillation, on each circuit of x2, x3 phase space there is an increase in amplitude in the interval [0, 2] and a decrease in the amplitude in the interval ["2, 0]. The actual value of these changes depends on the precise decimal values of the amplitudes. This enables us to compute analytically the expected amplitude distribution using the Irwin-Hall distribution, Fig. 7. V. DISSECTION OF THE SPACE OF INITIAL CONDITIONS Since two different paths are possible at each corner, there are, in principle, 24 or 16 different circuits around the square. In practice, not all are possible, as shown below. Based on the analysis above, we are able to compute the path that will be taken for any initial condition. Initial points from which one of the 16 circuits is taken form regions whose boundaries can also be computed. To illustrate this, we calculate the regions in "2 < x2 < 0, "9 < x3 < "7 at the point where x1 ¼ 1 and x4 ¼ 0 from which trajectories follow each of the possible types of path at each corner of a circuit. FIG. 7. Distribution of amplitudes starting from an initial amplitude of 5p on the edge between corner IV and corner I for 10 000 initial points for 30 circuits compared with the theoretically computed Irwin-Hall density for l ¼ 5p; r2 ¼ n=3, n ¼ 50. This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-8 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) FIG. 8. Dissection of the space of initial conditions. Labels (read right-toleft) refer to the type of path made by a single circuit through the four corners, beginning from the point (x1,x2,x3,x4) ¼ (1,x2,x3,0), where x2 and x3 correspond to a point in the graph. ABAA, for example, indicates that from the initial point, the trajectory takes path B only at corner III. Points in "1 < x2 < 0 take path A through Corner I, and thus, the coordinates correspond to the input points to Corner I, I;i xI;i 2 and x3 , as defined previously. However, points in "2 < x2 < "1 take path B through Corner I, and the input point occurs one time unit later, at x1 ¼ 0, x4 ¼ 1, and the coordinates of the initial point are thus x2 ¼ xI;i 2 " 1 and þ 1 (see Fig. 8). x3 ¼ xI;i 3 From Table VI, we can calculate the conditions on initial points that determine whether a trajectory passes through each Corner on path A or B. The conditions are listed in Table VIII and determine the regions shown in Fig. 8. In each case, if the condition for path A to be taken at a corner fails, then path B is taken. At integer values of each expression, i.e., boundary points between even and odd values of b2x2 þ x3 c, etc., path A and path B are equivalent in that the maps for the two paths coincide there. Finally, if path A is taken at Corner I, then the conditions in Table VIII I;i apply to xI;i 2 and x3 , but if path B is taken at Corner I, then I;i the conditions are all inverted for xI;i 2 and x3 . For example, if an initial point is in "2 < x2 < "1, then b2x2 þ x3 c I;i I;i I;i mod 2 ¼ 1 " b2xI;i 2 þ x3 cmod 2 ¼ d2x2 þ x3 emod 2. We now show that the set of points for which AAxx circuits occur indefinitely is of measure zero. We begin by analyzing the maps of initial regions that follow such a circuit twice. Fig. 8 shows (shaded) the subregions of the AAxx regions which map into an AAxx region. Thus, these initial points follow an AAxx circuit twice successively (not necessarily the same circuit: AAAA can be followed by AAAB). In these cases, either AAABAAAA or AAABAAAB, the composition of the maps for the two circuits give M & x2 x3 ' ¼ & 5 "4 4 "3 '& ' & ' x2 a þ ; "a x3 (4) where a is an integer depending on which of the shaded regions one is in. Consider, for example, the shaded region in Fig. 8 with vertices ("1.5, "7), ("1, "23/3), ("1, "7.5), ("1.4, "7). It lies in an AAAB region whose map is & ' & '& ' & ' x "3 "2 x2 "20 þ M2 2 ¼ : x3 x3 2 1 2 Its image is a polygonal region with vertices ("1.5, "8), ("5/3, "23/3), ("2, "7.5), ("1.8, "7.8), which lies in another AAAB region, whose map is & ' & '& ' & ' x "3 "2 x2 "22 þ : M1 2 ¼ x3 x3 2 1 4 Thus, the combined map, M ¼ M1M2, is (4) with a ¼ 34. TABLE VIII. Conditions on initial points for path A at each corner. If the condition fails, path B is taken. Corner I II III IV (if III took A) IV (if III took B) Condition for path A bx2 cmod 2 ¼ 1 b2x2 þ x3 cmod 2 ¼ 0 b3x2 þ 2x3 cmod 2 ¼ 1 b2x2 þ x3 cmod 2 ¼ 0 b10x2 þ 7x3 cmod 2 ¼ 1 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-9 Shahrear, Glass, and Edwards The map (4) has a line segment of fixed points given by x2 þ x3 ¼ "a/4 but restricted, of course, to the domain of definition of the map, i.e., the subregion of the particular AAAx region under consideration. In the case of AAAA regions, the line segment does not intersect the domain, so there is no line of fixed points. However, for AAAB regions, it does. For example, in the shaded region considered above, this line is x2 þ x3 ¼ "8.5, which—restricted to the shaded part of the AAAB region—is the line segment between ("1.5, "7) and ("1, "7.5). The matrix in (4) has eigenvalue 1 with algebraic multiplicity 2, but geometric multiplicity 1, and eigenvector ("1, 1). Thus, the one-dimensional eigenspace of any of the fixed points is simply the line containing the invariant segment. Points off this line are not fixed. However, amplitude remains unchanged under this doublecircuit map, since, letting (y2, y3) ¼ M(x2, x3), jy2 j þ jy3 j ¼ "y2 " y3 ¼ "5x2 " 4x3 " a þ 4x2 þ 3x3 þ a ¼ "x2 " x3 ¼ jx2 j þ jx3 j; using the fact that all coordinates are negative. Thus, points off the invariant line stay on their own line of slope of "1 but move steadily along it until they leave the shaded region on which the map is defined. The rate at which they move is constant but depends on distance from the invariant line. Thus, if x2 þ x3 ¼ "a=4 þ e, then y2 " x2 ¼ 5x2 þ 4x3 þa " x2 ¼ 4ð"a=4 þ eÞ þ a ¼ 4e. Thus, points off the invariant line move by a constant amount after each application of the double-circuit map until they leave the AAAB region and move into a region that is not of type AAxx. And, we have shown that the set of points that follow AAxx circuits indefinitiely exist but lie only on particular invariant line segments, which are of measure zero in the set of initial conditions, which is two-dimensional. The shaded regions in Fig. 8 amount to a proportion of 1/12 of the total area in the figure, and this proportion thus applies to the entire set of initial points (keeping away from the small amplitude region). Numerical simulations show that amplitude remains constant following two circuits about 7% of the time, which is a little less than 1/12 ) 8.333%. Why this discrepancy exists is not clear, though we have not shown that the density over the set of initial points remains uniform, so the proportion of area may not correspond to proportion of circuits on randomly chosen trajectories. While any initial point in the two-dimensional shaded regions return to the same amplitude after two circuits, there are also line segments (one-dimensional sets) from which trajectories return to the same amplitude after one circuit for some of the circuit types. An example is the segment from ("1, "8.4) to ("4/7, "9) in Fig. 8, which follows the circuit ABBA. But, these are sets of measure zero, and even there, the points are not fixed under the circuit map (except in the case AAAB discussed above) and thus leave the line segment after a few steps. VI. DISCUSSION This work investigated properties of a four dimensional ordinary differential equation based on a class of Chaos 25, 033103 (2015) equations that had previously been investigated as models of genetic networks, for example, see Refs. 2, 6,7, and 14. The current four dimensional example was discovered several years ago in random searches of differential equations in this class with four variables. The weird phase plane structures evident in Figure 2, led to the current investigation. A notable feature of this equation is that two variables form a negative feedback loop in which there is a decaying amplitude. In order to make the problem more tractable, we modified the original equations to generate a piecewise linear equation in which computing the dynamics reduces to an arithmetic problem. The resulting equations can be thought of as piecewise system of two variables that is subjected to a periodic forcing. In this respect, the current problem has a connection to classic studies of periodically forced oscillating systems such as the periodically forced van der Pol equation,15 except that the forced oscillator here is neutrally stable rather than a limit cycle. The modifications enabled us to develop a sharp analysis of the dynamics for high amplitude oscillations. A large number of related problems could be investigated. Notice that in Table III, there is an entry 0 in the top 4 rows designated x2(t þ 1). This 0 could have appeared in any of the top 3 rows and maintained the same qualitive feature that if x3 is low, then x2 will become active most of the time. Consequently, the current truth tables are but one representative of a class of 64 (8 ( 8) networks in which an oscillation of two variables drives a negative feedback loop of the other two variables. Moreover, since there are vast number of stable oscillations in the class of piecewise linear equations motivating the current study,6 the number of possible ways these can be hooked together is likewise vast. Would each possibility represent a new case, or is there some chance of a more general theory? A striking feature of the current system is the diffusive nature of the fluctuations of amplitude in a deterministic dynamical system. Our analysis here is similar to the analyses of diffusion in other deterministic systems in which it has been possible to compute analytically the probability distribution of the displacement.14,16–19 Our analysis of the reduced system shows that in general, trajectories cycle about the origin in the (x2, x3) plane, with amplitude diffusing like a random walk. In quantifying the diffusion, we found that although on each circuit, there is an upwards jump in amplitude in [0, 2] and a downwards jump in ["2, 0], the diffusion is not quite as fast as would be suggested by these jumps being uniformly random on those intervals. This is at least partially accounted for by the existence of regions of nonzero measure in which amplitude is unchanged after two circuits. Behaviour in the small amplitude region (amplitude % 2) can be understood at a microscopic level (see the Appendix), and a complete analysis would be of interest, though we have not done that here. In general, since amplitude is nonnegative, 0 acts as a reflecting barrier for the diffusion. Embedded in the (x2, x3) plane are special trajectories that are This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-10 Shahrear, Glass, and Edwards periodic (after one or two circuits, for example) and others that grow steadily in amplitude without bound, or shrink steadily to zero and then grow without bound. There is, in general, a sensitive dependence on initial conditions that can only be suspended for finite time intervals in the non-zero measure regions that return for a while to the same amplitude. Going back to the original network, from which the reduced system was derived, there are clearly some differences, despite the close approximation of the behaviour of the two systems in general. First, the reduced system is unbounded, while the original network is bounded because of the exponential decay term for each variable. Second, trajectories in the original network asymptotically collapse to the origin while fluctuating irregularly, but in the reduced system, they do not. This is a consequence of the shrinking amplitude of the x1, x4-cycle in the original network, as opposed to its constant amplitude in the reduced system. There are clear similarities, however. When the amplitude is small in comparison to the magnitude of the focal points (x2 and x3 coordinates), the piecewise-exponential time evolution of x2 and x3 in the original system is closely approximated by their piecewise-linear time evolution in the reduced system. Thus, as the amplitude shrinks in the original system, the reduced system becomes a better and better approximation, though the time unit defined by each leg of the x1, x4-cycle also shrinks, and thus, the scale of the reduced system shrinks. Numerically, the behavior of the two systems over finite time intervals appears very similar, and we believe that the analysis of sensitive dependence and diffusion in the reduced system provides a plausible explanation for their existence in the original network. This is an interesting example of deterministic diffusion that we have shown to be amenable to analysis, by means of a similar, reduced system on which calculations can be made explicitly. Although the current work was motivated initially by studies of genetic switching networks, it would be unrealistic to think that the current dynamics could be generated by real biological systems. The large derivatives make this implausible. On the other hand, building asynchronous switching systems with this logical structure is possible.20,21 Investigation of the dynamics in analogous systems with similar structures may lead to new classes of physical systems capable of generating flucutating dynamics. Chaos 25, 033103 (2015) APPENDIX A: TABLES OF INPUT OUTPUT COORDINATES The following tables give the output coordinates as a function of the input coordinates for each of the corners (in C;i which we take a ¼ jxC;i 2 j > 0 and b ¼ jx3 j > 0). We use these to construct the input-output maps for each corner, as well as the map from the output of one corner to the input of the next. To facilitate checking of the computations, numerical values from simulations are also given, illustrating a particular transition for each case. I;i Corner I: "1 < xI;i 2 < 0 and x3 < 0. The transition I;i I;i sequence for x1 ¼ 1; x4 ¼ 0 is given in Table IX. The tranI;i sition sequence for xI;i 1 ¼ 0; x4 ¼ 1 is given in Table X. The change in amplitude at Corner I is given by ðð1 " aÞ þ ð"1 þ 2a þ bÞÞ " ða þ bÞ ¼ 0. In vector notation, the input-output map for Corner I is given by xI;o ¼ CI xI;i þ DI ; $ % $ % 1 0 1 I I and D ¼ . where C ¼ 2 1 1 The map from the output of Corner I to the input to Corner II is given by I;o I;o xII;i 2 ¼ x2 " dx3 e; I;o I;o xII;i 3 ¼ x3 " dx3 e : II;i Corner II: xII;i 2 > 1 and "1 < x3 < 0. The transition II;i II;i sequence for x1 ¼ 0; x4 ¼ 1 is given in Table XI. The tranII;i sition sequence for xII;i 1 ¼ 1; x4 ¼ 0 is given in Table XII. The change in amplitude at Corner II is given by ðða þ 2b " 1Þ þ ð1 " bÞÞ " ða þ bÞ ¼ 0. In vector notation, the input-output map for Corner II is given by II II;i II xII;o % ¼ C x þ$D ; % 1 "2 "1 where CII ¼ and DII ¼ . 0 1 1 The map from the output of Corner II to the input to Corner III is given by $ II;o ¼ xII;o xIII;i 2 2 " bx2 c; II;o ¼ xII;o xIII;i 3 3 þ bx2 c : ACKNOWLEDGMENTS L.G. and R.E. thank NSERC for financial support. P.S. thanks the University of Bari Aldo Moro for financial support. APPENDIX In the Appendixes, we give several additional results. To facilitate checking of the computations, in Appendix A, we give tables of input and output coordinates for each corner of the x2, x3 phase plane. In Appendix B, we prove that the flow is expanding in the large amplitude case. Appendix C discusses special cases of the flows for small amplitudes and for integer values of the variables on the threshold planes. I;i TABLE IX. Corner I transitions-path A ða ¼ jxI;i 2 j; b ¼ jx3 jÞ. (a) x1 1 1"a 0 (b) Time 31939.50 31940.3638 31940.50 x2 x3 x4 "a 0 1"a "b "a " b 1 " 2a " b 0 a 1 x1 x2 x3 x4 1 0.1362 0 "0.8638 0 0.1362 "11.4381 "12.3018 "12.1656 0 0.8638 1 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-11 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) I;i TABLE X. Corner I transitions-path B ða ¼ jxI;i 2 j; b ¼ jx3 jÞ. (a) x1 0 "a "1 "a 0 a 1 (b) Time 31732.50 31733.2753 31733.50 31733.7246 31734.50 31735.27537 31735.50 TABLE XII. Corner II transitions-path B ða ¼ jx2II;i j; b ¼ jx3II;i jÞ. (a) x1 x2 x3 x4 "a 0 1"a 0 "a 0 1"a "b "a " b "1 " b "a " b "b "a " b 1 " 2a " b 1 1-a 0 a"1 "1 a"1 0 1 1"b 0 b"1 "1 b"1 0 (b) Time x1 x2 x3 x4 0 "0.775367 "1 "0.775367 0 0.775367 1 "0.775367 0 0.2246 0 "0.775367 0 0.2246 "11.27466 "12.0500 "12.27466 "12.0500 "11.27466 "12.0500 "11.8254 1 0.2246 0 "0.2246 "1 "0.2246 0 Corner III: 0 < xIII;i < 1 and xIII;i > 1. The transition 2 3 III;i III;i sequence for x1 ¼ "1; x4 ¼ 0 is given in Table XIII. The ¼ 0; xIII;i ¼ "1 is given in Table transition sequence for xIII;i 1 4 XIV. Table XIII suggests that the change in amplitude at Corner III on path A is ðð1 " aÞ þ ð1 þ bÞÞ " ða þ bÞ ¼ 2 " 2a. However, since this represents the high point on the squarewave in Figure 2, and we define the amplitude to be the lowest value of 3 consecutive points, the change in amplitude at Corner III on path A is actually "2a. The change in amplitude on path B is ðð1 " aÞ þ ð4a þ b " 3ÞÞ "ða þ bÞ ¼ "2 þ 2a. Thus, at Corner III, there is a decrease in amplitude that would be a random number between "2 and 0, provided a is a random number between 0 and 1. In vector notation, the input-output map for Corner III on path A is given by 32063.50 32063.8608 32064.500 32065.1392 32065.50 32065.8608 32066.500 x2 x3 x4 a aþb 1þa aþb a þ 2b aþb a þ 2b " 1 "b 0 1"b 0 "b 0 1"b 0 b 1 b 0 "b "1 x1 x2 x3 x4 1 0.6392 0 0.6392 "1 "0.6392 0 11.4675 11.8283 12.4675 11.82823 12.1891 11.8283 11.1891 "0.3608 0 0.6392 0 "0.3608 0 0.6392 0 0.3608 1 0.3608 0 "0.3608 "1 III;i xIII;o ¼ CIII þ DIII B x B ; $ % $ % 1 0 "1 III where CIII ¼ and D ¼ . B B 4 1 "3 The map from the output of Corner III to the input to Corner IV is given by ¼ xIII;o " bxIII;o xIV;i;AA 2 2 3 c; ¼ xIII;o " bxIII;o xIV;i;BA 2 2 3 c þ 2; ¼ xIII;o " bxIII;o xIV;i;AB 2 2 3 c " 2; ¼ xIII;o " bxIII;o xIV;i;BB 2 2 3 c; xIV;i ¼ xIII;o " bxIII;o 3 3 3 c; III;i xIII;o ¼ CIII þ DIII A x A ; $ % $ % 1 0 "1 III where CIII ¼ and D ¼ , while on path B, A A 0 1 1 it is given by where AB in the superscript, for example, indicates that path B is taken at Corner III and path A is taken at Corner IV. < "1 and 0 < xIV;i < 1. The transition Corner IV: xIV;i 2 3 IV;i sequence for x1 ¼ 1; xIV;i ¼ 0 is given in Table XV. The 4 ¼ 0; xIV;i ¼ 1 is given in Table transition sequence for xIV;i 1 4 XVI. Table XV suggests that the change in amplitude at TABLE XI. Corner II transitions-path A ða ¼ jx2II;i j; b ¼ jxII;i 3 jÞ. TABLE XIII. Corner III transitions-path A ða ¼ jx2III;i j; b ¼ jx3III;i jÞ. (a) x1 0 b 1 (b) Time 31758.50 31759.3254 31759.50 x2 x3 x4 (a) x1 x2 x3 x4 a aþb a þ 2b " 1 "b 0 1"b "1 b"1 0 "1 a"1 0 a 0 a"1 b aþb 1þb 0 "a "1 (b) Time x1 x2 x3 x4 0 0.8254 1 11.2246 12.0500 11.8754 "0.8254 0 0.1746 "1 "0.1746 0 32089.50 32089.6891 32090.500 x1 x2 x3 x4 "1 "0.81095 0 0.1891 0 "0.81095 11.6392 11.8283 12.6392 0 "0.18917 "1 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-12 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) TABLE XIV. Corner III transitions-path B ða ¼ jx2III;i j; b ¼ jx3III;i jÞ. IV;i TABLE XVI. Corner IV transitions-path B ða ¼ jxIV;i 2 j; b ¼ jx3 jÞ. (a) x1 x2 x3 x4 (a) x1 0 a 1 a 0 "a "1 a 0 a"1 0 a 0 a"1 b aþb 2a þ b " 1 3a þ b " 2 4a þ b " 2 3a þ b " 2 4a þ b " 3 "1 a"1 0 1"a 1 1"a 0 0 "b "1 "b 0 b 1 (b) Time (b) Time 42030.50 42030.8307 42031.50 0 42032.1693 42032.50 42032.8307 42033.500 x1 x2 x3 x4 0 0.3307 1 0.3307 0 "0.3307 "1 0.3307 0 "0.6693 0 0.3307 0 "0.6693 13.04336 13.374 12.7047 12.0354 12.3660 12.0354 11.3660 "1 "0.6693 0 0.6693 1 0.6693 0 Corner IV on path A is ða " 1 þ 1 " bÞ " ða þ bÞ ¼ "2b. However, since the entry point to corner IV in path A is at a high point on the square tooth wave in Figure 2, the change in amplitude is actually 2–2b. For the transition in Table XVI, the change in amplitude at Corner IV on path B is ð4b þ a " 1 þ 1 " bÞ " ða þ bÞ ¼ 2b. Thus, at corner IV, there is an increase in amplitude that is a random number between 0 and 2 if b is a random number between 0 and 1. In vector notation, the input-output map for Corner IV on path A is given by IV;i xIV;o ¼ CIV þ DIV A x A ; $ % $ % 1 0 1 IV IV where CA ¼ and DA ¼ , while on path B, 0 1 "1 it is given by IV;i xIV;o ¼ CIV þ DIV B x B ; $ % $ % 1 "4 1 IV where CIV ¼ and D ¼ . B B 0 1 "1 The map from the output of Corner IV to the input to Corner I is given by 42056.50 42056.8661 42057.50 42058.1339 42058.50 42058.8661 42059.50 0 (a) x1 1 1"b 0 (b) Time 31915.50 31916.06195 31916.50 x3 x4 "a b"a 1"a b 0 b"1 0 b 1 x1 x2 x3 x4 1 0.4381 0 "12.8638 "12.30184 "11.8638 0.56195 0 "0.4381 0 0.56195 1 x4 "a "a " b 1 " a " 2b "a " b "a " 2b "a " 3b 1 " a " 4b b 0 b"1 0 b 0 b"1 1 1"b 0 b"1 "1 b"1 0 x1 x2 x3 x4 0 "0.3661 "1 "0.3661 0 0.3661 1 "11.6693 "12.0354 "11.4015 "12.0354 "12.4015 "12.7675 "12.1336 0.3661 0 "0.6339 0 0.3661 0 "0.6339 1 0.6339 0 "0.6339 "1 "0.6339 0 IV;o þ dxIV;o xI;i 3 ¼ x3 2 e: APPENDIX B: PROOF OF EXPANSION Here, we prove that the iteration of return maps in our modified system is expanding unless a trajectory consists entirely of AAxx circuits. Lemma 1. & ' Let M be the class of 2 ( 2 real matrices 1 b M¼a satisfying the following properties: c d (1) (2) (3) (4) detðMÞ ¼ 1, jaj > 1, 1 , 0 < b < 1 " jaj 1 c > "1 þ jaj. M is closed under matrix multiplication. Proof: First, we observe that Let Then, x2 x3 IV;o xI;i " dxIV;o 2 ¼ x2 2 e; detð MÞ ¼ 1 M1 ¼ a1 IV;i TABLE XV. Corner IV transitions-path A ða ¼ jxIV;i 2 j; b ¼ jx3 jÞ. x2 & 1 c1 ' b1 ; d1 ) and M ¼ M1 M2 ¼ a with 1 : a2 & ' 1 b2 M2 ¼ a2 : c2 d2 d ¼ bc þ & 1 c ' b ; d a ¼ a1 a2 ð1 þ b1 c2 Þ; b1 ; a22 ð1 þ c2 b1 c2 Þ c ¼ c1 þ 2 : a1 ð1 þ b1 c2 Þ b ¼ b2 þ Note that (1 þ b1c2) > 0, since either c2 $ 0, in which case it is obvious, or c2 < 0 and jc2j < 1 and jb1j < 1. We check the four properties required to be in class M: This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-13 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) (1) detM ¼ detðM1 M2 Þ ¼ detðM1 ÞdetðM2 Þ ¼ 1 : ( ( )( )) (2) jaj¼ja1 jja2 jð1þb1 c2 Þ>ja1 jja2 j 1þ 1" ja11 j "1þ ja12 j ¼ja1 jþja2 j"1>1: (3) b > 0 since b2 > 0; b1 > 0, and ð1 þ b1 c2 Þ > 0 : Also, b1 þ b1 c2 Þ 1 b1 b1 ð1 " c2 ja2 jÞ " ja2 j <1" þ 2 ¼1þ ja2 j a2 ð1 þ b1 c2 Þ a22 ð1 þ b1 c2 Þ ja2 jðb1 " 1Þ ; <1þ 2 a2 ð1 þ b1 c2 Þ 1 since c2 > "1 þ ) 1 " c2 ja2 j < ja2 j: ja2 j b ¼ b2 þ a22 ð1 Thus $ % 1 "1 ja2 j 1 " ja1 j ; b<1þ a22 ð1 þ b1 c2 Þ since b1 < 1 " 1 ; ja1 j so b<1" 1 1 ¼1" : ja1 jja2 jð1 þ b1 c2 Þ jaj (4) c2 1 c2 þ > "1 þ ja1 j a21 ð1 þ b1 c2 Þ a21 ð1 þ b1 c2 Þ c2 ð1 þ ja1 jb1 Þ þ ja1 j : ¼ "1 þ a21 ð1 þ b1 c2 Þ c ¼ c1 þ Now, if c2 % 0, then b1 < 1 " ja11 j ) 1 þ ja1 jb1 < ja1 j, so c > "1 þ ðc2 þ 1Þja1 j 1 1 > "1 þ ¼ "1 þ ; 2 ð Þ ja1 jja2 j 1 þ b1 c2 jaj a1 ð1 þ b1 c2 Þ and if c2 > 0, then c2 ð1 þ b1 ja1 jÞ > 0, so ja1 j 1 and c > "1 þ again. c2 ð1 þ b1 ja1 jÞ þ ja1 j > ja1 j > ja j jaj 2 Thus, M 2 M: ! " be the closure of M and @M ¼ MnM " Let M be the boundary of M. Then, M 2 @M if and only if at least one of the 1 ; following equalities is satisfied: jaj ¼ 1; b ¼ 0; b ¼ 1 " jaj 1 c ¼ "1 þ jaj. " is closed under matrix multiplication. Lemma 2. M Proof: This can be shown by an obvious modification to the proof of Lemma 1 or by continuity. ! Now the following result can be checked by looking at boundary cases in the proof of Lemma 1. Lemma 3. If M1 2 M and M2 2 @M then M1 M2 2 M and M2 M1 2 M. A matrix, M, has an expanding direction if its spectral radius, qðMÞ, is greater than 1. Lemma 4. M 2 M ) qðMÞ > 1. " ) detðMÞ ¼ 1, which implies that both Proof: M 2 M eigenvalues are real. We can label the eigenvalues such that qðMÞ ¼ jk1 j $ jk2 j. Then, k2 ¼ k11 and jtraceðMÞj ¼ jk1 þ k11 j ¼ jk1 j þ jk11 j. Thus, qðMÞ ¼ 1 if and only if jtraceðMÞj ¼ 2, and otherwise, qðMÞ > 1 and jtraceðMÞj > 2. Now, M 2 M ) $ %* * * * * 1 ** ** 1** * ð Þ ð Þ jtrace M j ¼ *a þ a bc þ 2 * ¼ *a 1 þ bc þ * a a 1 ¼ jajð1 þ bcÞ þ ; jaj since ð1 þ bcÞ > 0 : Thus, $ $ %$ %% 1 1 1 jtraceð MÞj > jaj 1 þ 1 " "1 þ þ jaj jaj jaj ! 2 1 1 þ " ¼ 2: ¼ jaj jaj jaj2 jaj Thus, qðMÞ > 1: ! Denote the matrices in Table VII in the natural way as MAA ; MBA ; MAB , and MBB. Then MAA 2 @M, while MBA, k 2 @M for MAB, and MBB 2 M. One can check that MAA k ¼ 1,2,3,…. However, Lemma 3 implies that if M is any other multiple product (word) of the above four matrices, k Þ ¼ 2 for then M 2 M. One can check that traceðMAA k k ¼ 1,2,3,…, so qðMAA Þ ¼ 1, but if M is any other word made from the above four matrices, then by Lemma 4, qðMÞ > 1, and it has an expanding direction. APPENDIX C: SMALL AMPLITUDE BEHAVIOR AND INTEGER POINTS When the amplitude of a trajectory is large, many transitions of x1 and x4 occur between transitions involving x2 or x3 (the paths between the corners). When the amplitude is small, however, it is possible for multiple transitions of x2 and x3 to occur between successive transitions of x1 or x4, i.e., between clock steps of our process. Thus, it is helpful to consider the behavior in the x2, x3-plane during each of the four states of (X1, X4). This is illustrated in Fig. 9, for X1, X4 ¼ 00, 10, 11, and 01 in the four panels. For example, over the unit time interval when X1 ¼ 1 and X4 ¼ 0 (i.e., x1 > 0, x4 < 0) the unit-length trajectory segment in (x2, x3) follows a constant-amplitude cycle counterclockwise around the origin. Far from the origin in one quadrant, this is just a unit-length line segment without transitions, but close to the origin, it may involve many transitions and even many cycles. Note that at the origin, (x2, x3) ¼ (0, 0), during this phase of the x1, x4-cycle, the flow is not well defined, but the only consistent solution is one that remains constant at the origin. In the subsequent phase, X1, X4 ¼ 11, of course, the trajectory then leaves the origin, going to (x2, x3) ¼ (1, 1) over the next time step. The flows shown in Fig. 9 make it clear that there is no non-uniqueness of trajectories, even if x2 or x3 is 0 at a clock step (i.e., an integer time point when x1 or x4 is 0—boundaries between A and B paths around corners in the large amplitude This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-14 Shahrear, Glass, and Edwards Chaos 25, 033103 (2015) FIG. 9. Unit length trajectories in each of the four phases of the (X1, X4) cycle. (a) (X1,X4) ¼ (0,0), (b) (X1,X4) ¼ (1,0), (c) (X1,X4) ¼ (1,1), (d) (X1,X4) ¼ (0,1). For example, (X1,X4) ¼ (0,0) for (x1, x4) between ("1, 0) and (0, "1). case are such points, for example). The region around the origin in the x2, x3-plane can be partitioned into regions in which different behaviors (sequences of transitions) occur. The flow can be complicated, and if we define ‘small amplitude’ to be amplitude % 2 (keeping in mind our definition of amplitude as the minimum of jx2j þ jx3j over three succesive time steps), it is possible for some trajectories to remain in the small amplitude region for a long time, though most are eventually ejected again. In particular, a trajectory starting from (x1, x4) ¼ (1, 0) (as in Fig. 8) but (x2, x3) ¼ (0, "1) is periodic and has amplitude 1 indefinitely, though it makes several brief excursions to points where jx2j þ jx3j ¼ 3. In fact, it can be shown that starting at (x2, x3) points with integer coordinates and odd amplitude are always periodic after a single circuit, even in the small amplitude case. The point ("1, "8) in Fig. 8 is an example, as are the four corner points in that figure, ("2, "9), ("2, "7), (0, "7), and (0, "9). Furthermore, integer points with even amplitude return after one circuit with amplitude increased or decreased by 2, depending on whether x2 (or x3) is individually even or odd at a particular phase. For example, in Fig. 8, the point ("2, "8) returns to ("2, "10), while the point ("1, "9) returns to ("1, "7). Thus, starting at the phase (x1,x4) ¼ (1,0), points (x2,x3) ¼ (2j,2k) have trajectories that steadily grow in amplitude (by 2 each circuit), indefinitely, while points (2j þ 1,2k þ 1) have trajectories that steadily decrease in amplitude (by 2 each circuit), until they reach amplitude 0, where they stay for one time step (while X1,X4 ¼ 10) and thus shift phase, so that at (x1,x4) ¼ (1,0), they now have x2 and x3 even, so the amplitude then grows without bound. Integer points clearly do not display diffusive behaviour. The line segments that are fixed under two circuits of the map, discovered above, are also clearly not diffusive (these are points with half-integer amplitude). But these trajectories are special, and arise from measure zero sets of initial conditions. With regard to the diffusion in amplitude, we found in Section IV, there is clearly a lower bound at 0, and this acts as a reflecting barrier in general. Amplitude can diffuse down towards 0, (small amplitude being defined to be % 2), but will then generally diffuse up again, though it may also remain near 0 for long periods of time. 1 S. A. Kauffman, The Origins of Order: Self-organization and Selection in Evolution (New York, Oxford University Press, 1993). 2 L. Glass and S. Kauffman, “The logical analysis of continuous, non-linear biochemical control networks,” J. Theor. Biol. 39, 103–29 (1973). 3 L. Glass and J. S. Pasternack, “Prediction of limit cycles in mathematical models of biological oscillations,” Bull. Math. Biol. 40, 27–44 (1978). 4 T. Mestl, E. Plahte, and S. Omholt, “Periodic solutions in systems of piecewise-linear differential equations,” Dyn. Syst. 10, 179–193 (1995). 5 E. Plahte and S. Kjøglum, “Analysis and generic properties of gene regulatory networks with graded response functions,” Physica D 201, 150–176 (2005). 6 R. Wilds and L. Glass, “An atlas of robust, stable, high-dimensional limit cycles,” Int. J. Bifurc. Chaos 19, 4055–4096 (2009). 7 T. Mestl, C. Lemay, and L. Glass, “Chaos in high-dimensional neural and gene networks,” Physica D 98, 33–52 (1996). This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 033103-15 8 Shahrear, Glass, and Edwards R. Edwards, “Chaos in neural and gene networks with hard switching,” Differ. Equations Dyn. Syst. 9, 187–220 (2001). 9 L. Lu and R. Edwards, “Structural principles for complex dynamics in glass networks,” Int. J. Bifurc. Chaos 21, 237–254 (2011). 10 P. Shahrear, L. Glass, R. Wilds, and R. Edwards, “Dynamics in piecewise linear and continuous models of complex switching networks,” Math. Comput. Simul. (in press). 11 P. Hall, “The distribution of means for samples of size n drawn from a population in which the variate takes values between 0 and 1, all such values being equally probable,” Biometrika 19, 240–245 (1927). 12 J. O. Irwin, “On the frequency distribution of the means of samples from a population having any law of frequency with finite moments, with special reference to Pearson’s Type II,” Biometrika 19, 225–239 (1927). 13 L. Glass and J. S. Pasternack, “Stable oscillations in mathematical models of biological control systems,” J. Math. Biol. 6, 207–223 (1978). Chaos 25, 033103 (2015) 14 T. Mestl, R. Bagley, and L. Glass, “Common chaos in arbitrarily complex feedback networks,” Phys. Rev. Lett. 79, 653 (1997). M. Levi, Qualitative Analysis of the Periodically Forced Relaxation Oscillations (Memoirs of the American Mathematical Society, 1981), p. 244. 16 T. Geisel and J. Nierwetberg, “Onset of diffusion and universal scaling in chaotic systems,” Phys. Rev. Lett. 48, 7 (1982). 17 R. Klages and J. Dorfman, “Simple deterministic dynamical systems with fractal diffusion coefficients,” Phys. Rev. E 59, 5361 (1999). 18 R. Klages, Microscopic Chaos, Fractals and Transport in Nonequilibrium Statistical Mechanics (World Scientific, Singapore, 2007). 19 J. Lei and M. C. Mackey, “Deterministic Brownian motion generated from differential delay equations,” Phys. Rev. E 84, 041105 (2011). 20 J. Mason, P. S. Linsay, J. J. Collins, and L. Glass, “Evolving complex dynamics in electronic models of genetic networks,” Chaos 14, 707–715 (2004). 21 D. P. Rosin, D. Rontani, D. J. Gauthier, and E. Sch€ oll, “Experiments on autonomous Boolean networks,” Chaos 23, 025102 (2013). 15 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.204.177.148 On: Wed, 11 Mar 2015 11:32:36 View publication stats