A Proofs

Shaping Pulses to Control Bistable Biological Systems.

Abstract

In this paper we study how to shape temporal pulses to switch a bistable system between its stable steady states. Our motivation for pulse-based control comes from applications in synthetic biology, where it is generally difficult to implement real-time feedback control systems due to technical limitations in sensors and actuators. We show that for monotone bistable systems, the estimation of the set of all pulses that switch the system reduces to the computation of one non-increasing curve. We provide an efficient algorithm to compute this curve and illustrate the results with a genetic bistable system commonly used in synthetic biology. We also extend these results to models with parametric uncertainty and provide a number of examples and counterexamples that demonstrate the power and limitations of the current theory. In order to show the full potential of the framework, we consider the problem of inducing oscillations in a monotone biochemical system using a combination of temporal pulses and event-based control. Our results provide an insight into the dynamics of bistable systems under external inputs and open up numerous directions for future investigation.

\IEEEoverridecommandlockouts

1 Introduction

In this paper we investigate how to switch a bistable system between its two stable steady states using external input signals. Our main motivation for this problem comes from synthetic biology, which aims to engineer and control biological functions in living cells [3]. Most of current research in synthetic biology focusses on building biomolecular circuits inside cells through genetic engineering. Such circuits can control cellular functions and implement new ones, including cellular logic gates, cell-to-cell communication and light-responsive behaviours. These systems have enormous potential in diverse applications such as metabolic engineering, bioremediation, and even the energy sector [4].

Several recent works [5, 6, 7] have showcased how cells can be controlled externally via computer-based feedback and actuators such as chemical inducers or light stimuli [8, 9]. An important challenge in these approaches is the need for real-time measurements, which tend to be costly and difficult to implement with current technologies. In addition, because of technical limitations and the inherent nonlinearity of biochemical interactions, actuators are severely constrained in the type of input signals they can produce. As a consequence, the input signals generated by traditional feedback controllers (e.g. PID or model predictive control) may be hard to implement without a significant decrease in control performance.

In this paper we show how to switch a bistable system without the need for output measurements. We propose an open-loop control strategy based on a temporal pulse of suitable magnitude and duration :

(1)

Our goal is to characterise the set of all pairs that can switch the system between the stable steady states and the set of all pairs that cannot. We call these sets the switching sets and a boundary between these sets the switching separatrix. The pairs close to the switching separatix are especially important in synthetic biology applications, as a large or a large can trigger toxic effects that slow down cell growth or cause cell death.

In a previous paper [1], we showed that for monotone systems the switching separatrix is a monotone curve. This result was therein extended to a class of non-monotone systems whose vector fields can be bounded by vector fields of monotone systems. This idea ultimately leads to robustness guarantees under parametric uncertainty. These results are in the spirit of [10, 11, 12], where the authors considered the problem of computing reachability sets of a monotone system. Some parallels can be also drawn with [13, 14], where feedback controllers for monotone systems were proposed.

Contributions. In the present paper we provide the first complete proof of our preliminary results in [1] and extend them in several directions. We formulate necessary and sufficient conditions for the existence of the monotone switching separatrix for non-monotone systems. Although it is generally hard to use this result to establish monotonicity of the switching separatrix, we use it to prove the converse. For example, we show that for a bistable Lorenz system the switching separatrix is not monotone. We then generalise the main result of [1] by providing conditions for the switching separatrix to be a graph of a function. We also discuss the relation between bifurcations and the mechanism of pulse-based switching, which provides additional insights into the switching problem. We use this intuition to show and then explain the failure of pulse-based control on an HIV viral load control problem [15]. We proceed by providing a numerical algorithm to compute the switching separatrices for monotone systems. The algorithm can be efficiently distributed among several computational units and does not explicitly use the vector field of the model. We evaluate the theory and computational tools on the bistable LacI-TetR system, which is commonly referred to as a genetic toggle switch [16]. Genetic toggle switches are widely used in synthetic biology to trigger cellular functions in response to extracellular signals [3, 17].

We complement our theoretical findings with several observations that illustrate limitations of the current theory and highlight the need for deeper investigations of bistable systems. For example, we show that for a toxin-antitoxin system [18], the switching separatrix appears to be monotone, even though the system does not appear to be monotone. Finally, in order to demonstrate the full potential of pulse-based control, we consider the problem of inducing an oscillatory behaviour in a generalised repressilator system [19].

Organisation. The rest of the paper is organised as follows. In Section 2 we cover the basics of monotone systems theory, formulate the problem in Subsection 2.1, and provide an intuition into the mechanism of pulse-based switching for monotone systems in Subsection 2.2. We also provide some motivational examples for the development of our theoretical results. In Section 3 we formulate the theoretical results and in Section 4 we present the computational algorithm, which we evaluate in Section 5 on the LacI-TetR system. In Section 6, we provide counterexamples and an application of inducing oscillations in a generalised repressilator system. All the proofs are found in the Appendix.

Notation. Let stand for the Euclidean norm in , stand for a topological dual to , stand for the relative complement of in , stand for the interior of the set , and for its closure.

2 Preliminaries

Consider a single input control system

(2)

where , , , and belongs to the space of Lebesgue measurable functions with values from . We say that the system is unforced, if . We define the flow map , where is a solution to the system (2) with an initial condition and a control signal . We consider the control signals in the shape of a pulse, that is signals defined in (1) with the set of admissible and denoted as .

In order to avoid confusion, we reserve the notation for the vector field of non-monotone systems, while systems

(3)
(4)

denote so-called monotone systems throughout the paper. In short, monotone systems are those which preserve a partial order relation in initial conditions and input signals. A relation is called a partial order if it is reflexive (), transitive (, implies ), and antisymmetric (, implies ). We define a partial order through a cone as follows: if and only if . We write , if the relation does not hold. We will also write if and , and if . Similarly we can define a partial order on the space of signals : if for all . We write , if and for all . Partial orders also induce some geometric properties on sets. A set is called p-convex if for every , in such that , and every we have that .

Definition 1

The system (3) is called monotone on with respect to the partial orders , , if for all and such that and , we have for all . If additionally, , or implies that for all , then the system is called strongly monotone.

In general, it is hard to establish monotonicity of a system with respect to an order other than an order induced by an orthant (e.g., positive orthant ). Hence throughout the paper, by a monotone system we actually mean a monotone system with respect to a partial order induced by an orthant. A certificate for monotonicity with respect to an orthant is referred to as Kamke-Müller conditions [20].

Proposition 1 ([20])

Consider the system (3), where is differentiable in and and let the sets , be p-convex. Let the partial orders , be induced by , , respectively, where , for some , in . Then

if and only if the system (3) is monotone on with respect to ,

If we consider the orthants , , then the conditions above are equivalent to checking if for all such that for some , and all we have .

2.1 Problem Formulation

We confine the class of considered control systems by making the following assumptions:

  1. Let in (2) be continuous in on . Moreover, for each compact sets and , let there exist a constant such that for all and .

  2. Let the unforced system (2) have two stable steady states in , denoted as and ,

  3. Let , where stands for the domain of attraction of the steady state for of the unforced system (2),

  4. For any with finite and let belong to . Moreover, let the sets

    be non-empty.

Assumption A1 guarantees existence, uniqueness and continuity of solutions to (2), while Assumptions A2–A3 define a bistable system on a set controlled by pulses. In Assumption A4 we define the switching sets: the set , which contains all pairs that switch the system, and the set , which contains all pairs that do not. The boundary between these sets is called the switching separatrix. In the rest of the paper, we focus on the control problem of estimating the switching sets.

2.2 Mechanism of Pulse-Based Switching

The general problem of switching a bistable system with external inputs is amenable to an optimal control formulation. However, in applications such as synthetic biology, optimal control solutions can be very hard to implement due to technical limitations in actuators and output measurements. Additionally, the solution of this optimal control problem may be technically challenging. Hence applying open-loop pulses can be a reasonable solution, if we can guarantee some form of robustness. As we shall see later, our results show that for monotone systems, pulse-based switching is computationally tractable and robust towards parameter variations.

Before presenting our main results, we first provide an intuitive link between monotonicity and the ability to switch a system with temporal pulses. If we consider constant inputs and regard as a bifurcation parameter, we have the following result.

Proposition 2

Let the system (3) satisfy Assumptions A1–A4 and be monotone on with respect to , . Let be such that all pairs for , and any finite positive . Let also and . Then

  1. If then ,

  2. If then and .

  3. The function is discontinuous at .

Figure 1: A schematic depiction of the evolution of the stable nodes , and the saddle with respect to in the genetic toggle switch system. By slow manifold we mean a manifold connecting stable equillibria and a saddle. The arrows show the direction of the equillibria movements with increasing . At the equilibria and collide resulting in a saddle-node bifurcation preserving only .

The proof of the proposition is in the Appendix. In many applications, the functions , are simply evolutions of the steady states , with respect to the parameter , respectively. Hence, statement (1) of Proposition 2 shows how the steady states move with respect to changes in . Statement (2) ensures that there are at least two distinct asymptotically stable equilibria for . Finally, statement (3) indicates that the system undergoes a bifurcation for . The particular type of the bifurcation will depend on a specific model. Next we investigate further aspects of this result with some examples of monotone and non-monotone bistable systems.

Example 1: LacI-TetR Switch. The genetic toggle switch is composed of two mutually repressive genes LacI and TetR and was a pioneering genetic system for synthetic biology [16]. We consider its control-affine model, which is consistent with a toggle switch actuated by light induction [9]:

(5)

where the parameters have the following values

(6)

In the model (5), represents the concentration of each protein, whose mutual repression is modelled via a rational function. The parameters and represent the repression thresholds, whereas and model the basal synthesis rate of each protein. The parameters and are the degradation rate constants, , are called Hill (or cooperativity) parameters, and , describe the strength of mutual repression. By means of Proposition 1 we can readily check that the model is monotone on for all nonnegative parameter values. It can be verified by direct computation that the system satisfies Assumptions A1–A4 with . It can be also shown that the unforced system is strongly monotone in using the results in [21].

With the chosen parameter values, we numerically found a bifurcation to occur at . For the system has two stable nodes and a saddle. We observe that for all , and therefore we conclude that the system undergoes a saddle-node bifurcation, as illustrated in Figure 1.

Example 2: Lorenz system. Consider a system

with parameters , , , which is non-monotone and bistable with two stable foci. Numerical computation of the sets and in Figure 2 suggests that the switching separatrix is not monotone. We will revisit this conclusion in the next section using our theoretical results.

Figure 2: Switching sets for the Lorenz system. We simulated the Lorenz system for pairs taken from a mesh grid. The green and red crosses correspond to the pairs that switched or not switched the system, respectively.

Example 3: HIV viral load control problem. In [15] the authors considered the problem of switching from a “non-healthy” () to a “healthy” () steady state by means of two control inputs ( and ) that model different drug therapies. Due to space limitations we refer the reader to [15] for a description of the model. It can be shown that both steady states are stable foci and that the model is non-monotone. Although the system can be switched with non-pulse control signals [15], using extensive simulations we were unable to find a combination of pulses in and switching the system.

As in the case of monotone bistable system, we found a bifurcation with respect to constant control signals and . More specifically, we fixed , and numerically found a bifurcation at . The major difference between this case and the monotone system case (Example 1) is that the steady state lies the domain of attraction of . Hence if we stop applying the constant control signal we regress back to the initial point . Furthermore, with increasing the steady state is moving towards the origin, which also lies in the domain of attraction of . This makes pulse-based switching very difficult, if not impossible.

3 Theoretical Results

In [1] we showed that the switching separatrix of a monotone bistable system is non-increasing. Here we present a generalisation of this result by formulating necessary and sufficient conditions for the switching separatrix to be monotone, the proof of which is found in the Appendix.

Theorem 1

Let the system (2) satisfy Assumptions A1–A4. Then the following properties are equivalent:

  1. If belongs to for all , then belongs to for all , and for all , such that , .

  2. The set is simply connected. There exists a curve , which is a set of maximal elements of in the standard partial order. Moreover, the curve is such that for any and , for .

Theorem 1 shows that the computation of the set is reduced to the computation of a curve . This result also provides a connection between the geometry of domains of attraction of the unforced system and the switching separatrix. As shown next, Theorem 1 can also be used to establish non-monotonicity of the switching separatrix.

Remark 1 (Lorenz system revisited)

Consider the Lorenz system from the previous section and three different pulses with , , , and . Numerical solutions with increased accuracy show that the flows and converge to , whereas converges to . Application of Theorem 1 proves that the switching separatrix is not monotone.

The major bottleneck in the direct application of Theorem 1 is the verification of condition (1), which is generally computationally intractable. For example, condition (1) is satisfied if the partial order is preserved for control signals. That is for any , it should follow that for all . Although this property is weaker than monotonicity, it is not clear how to verify it. Monotonicity, on the other hand, is easy to check and implies condition (1) in Theorem 1. This is used in the following result.

Theorem 2

Let the system (3) satisfy Assumptions A1–A4 and be monotone on .

  1. The set is simply connected. There exists a curve , which is a set of maximal elements of in the standard partial order. Moreover, the curve is such that for any and , for .

  2. The set is simply connected. There exists a curve , which is a set of minimal elements of in the standard partial order. Moreover, the curve is such that for any and , for .

  3. Let the system (3) be strongly monotone and be the separatrix between the domains of attractions and of the unforced system (3). Let additionally be an unordered manifold, that is, there are no , in such that . Then for all and the curve is a graph of a monotonically decreasing function.

We state implicitly in Theorem 2, that if , then the flow does not converge to or , since it may end up on the separatrix . We note that our computational procedure presented in Section 4 does not require that or that , are graphs of functions. Hence we treat point (3) in Theorem 2 as a strictly theoretical result, but remark that sufficient conditions for the separatrix to be unordered are provided in Theorem 2.1 in [22]. The most relevant condition to our case is that the flow of the unforced system is strongly monotone, which we also assume in Theorem 2.

Besides , there are other pathological cases. For example, applying constant input control signals typically results in a system (2) with a different set of steady states than or . Moreover, the number of equilibria may be different. Hence, with the set typically does not contain the limiting control signal . If the set of pairs resulting in these pathological cases is not measure zero, then the sets and are not equal, which can complicate the computation of the switching sets. However, in many practical applications, the sets and appear to be equal. Therefore in order to simplify the presentation we study only the properties of the set .

If the system to be controlled is not monotone, then the curve may not be monotone, which is essential for our computational procedure. Instead, we estimate inner and outer bounds on the switching set provided that the vector field of the system can be bounded from above and below by vector fields of monotone systems. This is formally stated in the next result, while the proof is in the Appendix.

Theorem 3

Let systems (2), (3), (4) satisfy Assumptions A1–A4. Let , the systems (3) and (4) be monotone on and

(7)

Additionally assume that the stable steady states , , , satisfy

(8)
(9)

Then the following relations hold:

(10)

The technical conditions in (8), (9) are crucial to the proof and are generally easy to satisfy. An illustration of these conditions is provided in Figure 3. Checking the condition (9) reduces to the computation of the stable steady states, as does checking the condition (8). Indeed, to verify that belongs to the intersection of , , , we check if the trajectories of the systems (3), (4) initialised at with converge to and , respectively, which is done by numerical integration. The computation of stable steady states can be done using the methods from [23].

Figure 3: A schematic depiction of the conditions (8) and (9). The condition (8) ensures that all the steady states lie in the intersection of the corresponding domains of attractions (violet area). The steady state cannot lie in the dashed blue box due to condition (9).

In some applications, we need to find a subset of the pairs that switch the system (2) from to . Due to the inclusion , existence of the system (3) allows to do that. In this case, we are only interested in finding the system (3), hence the condition (9) is not required and the condition (8) is transformed to , .

Remark 2

The proofs of Theorems 2 and 3 are adapted in a straightforward manner to the case when systems are monotone with respect to an order induced by an arbitrary cone , and the order induced by . In examples, however, we always assume that is an orthant, since it is hard to check monotonicity of the system with respect to an arbitrary cone.

Theorem 3 also provides a way of estimating the switching set under parametric uncertainty, which is stated in the next corollary.

Corollary 1

Consider a family of systems with a vector of parameters taking values from a compact set . Let the systems satisfy Assumptions A1–A4 for every in . Assume there exist parameter values , in such that the systems and are monotone on , where and

(11)

for all . Let also

(12)
(13)

for all in . Then the following relation holds:

(14)

The proof follows by setting and and noting that the conditions in (12), (13) imply the conditions in (8), (9) in the premise of Theorem 3.

Theorem 3 states that if the bounding systems (3), (4) can be found, the switching sets , can be estimated, thereby providing approximations on the switching set . In what follows we provide a procedure to find monotone bounding systems if the system (2) is near-monotone, meaning that by removing some interactions between the states the system becomes monotone (see [24] for the discussion on near-monotone systems).

Let there exist a single interaction which is not compatible with monotonicity with respect to an order induced by , and let this interaction be between the states and . This happens if, for example, the -th entry in the Jacobian is smaller or equal to zero. A monotone system can be obtained by replacing the variable with a constant in the function , which removes the interaction between the states and . If the set is bounded then clearly we can find and such that for all . If the set is not bounded, then we need to estimate the bounds on the intersection of and the reachability set starting at for all admissible pulses. Let for all , , and . It is straightforward to show that , and are monotone systems and their vector fields are bounding the vector field from below and above, respectively. Note that in order to apply Theorem 3 we still need to check if these bounding systems satisfy Assumptions A1–A4.

In the case of Corollary 1, the procedure is quite similar. If the system is monotone for all parameter values , then we can find , if there exists a partial order in the parameter space. That is a relation such that for parameter values and satisfying we have that

If a partial order is found, the values and are computed as minimal and maximal elements of in the partial order . This idea is equivalent to treating parameters as inputs and showing that the system is monotone with respect to inputs and .

4 Computation of the Switching Separatrix

The theoretical results in Section 3 guarantee the existence of the switching separatrix for monotone systems, but in order to compute we resort to numerical algorithms.

Given a pair we can check if this pair is switching the system using simulations (that is, numerically integrating the corresponding differential equation). If the curve is a monotone function, then for every there exists a unique pulse magnitude . Let be such that for all . Clearly, for every we can compute the corresponding using bisection. We start the algorithm by computing the value corresponding to . Due to monotonicity of the switching separatrix, the minimal switching magnitude for the pulse length is smaller or equal to . Therefore, we can save some computational effort by setting the upper bound on the computation of equal to . The computation of the pairs can be parallelised by setting the same upper bound on , , , where is the number of independent computations. The procedure is summarised in Algorithm 1, where and are the sets of pairs approximating the switching separatrix from below and above, respectively.

1:Inputs: The system with initial state , final state , tolerance , simulation time , a grid , an upper bound on the magnitude , the number of used processors .
2:Outputs: finite sets and
3:for  do
4:      Set , ,
5:      for  do
6:            while  do
7:                 
8:                 Set
9:                 if  then
10:                       
11:                 else
12:                       
13:                 end if
14:            end while
15:      end for
16:      
17:      
18:end for
Algorithm 1 Bisection Algorithm for Computation of the Switching Separatrix
Figure 4: Illustration of the error of computation of the switching separatrix between the values , . The black curve is the switching separatrix to be computed, the red and green circles are the upper and lower bounding points, respectively. The switching separatrix should lie between the coloured regions due to its monotonicity. The values and are the largest height and width of boxes inscribed between the coloured regions, respectively.

In order to evaluate the error of computing the switching separatrix consider Figure 4. According to the definitions in the caption of Figure 4 we define the relative error of the approximation as

Note that, even if the green and red circles lie very close to each other the relative error can be substantial. In numerical simulations we use a logarithmic grid for , which yields a significantly lower relative error in comparison with an equidistant grid. This can be explained by an observation that in many numerical examples appears to be an exponentially decreasing curve.

1:Inputs: The system with initial state , final state , total number of samples , simulation time , lower and upper bounds on , and respectively, the numbers , , probability distribution
2:Outputs: sets and
3:Compute and using bisection for values and
4:Set
5:Set
6:for  do
7:      Compute the values , , and the corresponding boxes , .
8:      Generate samples in each of the boxes and using a probability distribution
9:      Generate randomly samples
10:      for  do
11:            Check if the samples switch the system
12:      end for
13:      Update and prune the sets ,
14:end for
Algorithm 2 Computation of Switching Separatrix Based on Random Sampling

There are a few drawbacks in Algorithm 1. Firstly, it requires a large number of samples. Secondly, the choice of the grid is not automatic, which implies that for switching separatrices with different geometry the relative error on the same grid may be drastically different. Finally, the algorithm relies on the assumption that is a graph of a monotone function, which may not be true. In order to overcome these difficulties, we have derived Algorithm 2 based on random sampling, which converges faster than Algorithm 1, has higher sample efficiency, does not require a predefined grid and the graph assumption. Some of the steps in Algorithm 2 require additional explanation:

Step 7. Find two boxes: the box with the maximal height (denoted as ) and the box with the maximal width (denoted as ) that can be inscribed between the coloured regions as depicted in Figure 4.

Step 9. Generate samples of using a probability distribution between and . For every generate a value using a distribution such that lies in the area between the coloured regions. Repeat this step by first generating between and using a distribution , and then generating for every generated in the area between the coloured regions.

Step 13. First, we update the sets , by adding the samples that do not switch and switch the system, respectively. Then if there exist two pairs and in the set (resp., ) such that and , then delete the pair from the set (resp., the pair from the set ).

Note that Step 11 is the most computationally expensive part of the algorithm and its computation is distributed into independent tasks. In our implementation, we chose as a Beta distribution with parameters and and adjusted the support to a specific interval. Note that the set between the coloured regions is getting smaller with every generated sample, hence the relative error of Algorithm 2 is a non-increasing function of the total number of samples. In fact, numerical experiments show that this function is on average exponentially decreasing. After the sets and are generated one can employ machine learning algorithms to build a closed form approximation of a switching separatrix (e.g., Sparse Bayesian Learning [25]; see also [26], [27] for efficient algorithms).

Figure 5: Average error against total number of generated samples. The curves corresponding to Algorithm 1 are computed by a single run of the algorithm. The curves corresponding to are averages over ten runs of Algorithm 2, while the curves for are the averages over twenty runs of Algorithm 2. Recall that for Algorithm 2.

5 Illustration of Theoretical Results on the LacI-TetR Switch

5.1 Evaluation of the Computational Algorithm

Here we will compare Algorithm 1 and 2 with different parameter values, as well as their distributed implementations on the LacI-TetR switch introduced in Subsection 2.2. Note that Algorithm 2 does not depend explicitly on the dynamics of the underlying system, but depends only on the generated pairs . Therefore, the convergence and sample efficiency results presented here will be valid for a broad class of systems. In Figure 5, we compare the error against the total number of generated samples. Since checking if a sample switches the system is the most expensive part of both algorithms, the total number of samples reflects the computational complexity. In the case of Algorithm 2 with the randomisation level is not high, hence an average over ten runs is sufficient to demonstrate the average behaviour of this algorithm. Note that both curves corresponding to Algorithm 2 with outperform the curves corresponding to Algorithm 1 in terms of accuracy.

Algorithm
Alg. 1 with
Alg. 1 with
Alg. 2 with ,
Alg. 2 with ,
Alg. 2 with ,
Alg. 2 with ,
Table 1: Sample efficiency in percent. In the notation , stands for the emperical mean, and for the emperical standard deviation. Recall that for Algorithm 2.

Some computational effort in Algorithm 2 goes into computing the error. However, this effort appears to be negligible in comparison with numerically solving a differential equation for a given pair even for such a small system as the toggle switch. We run the simulations on a computer equipped with Intel Core i7-4500U processor and 8GB of RAM. Using the centralised version of Algorithm 2 we achieved on average a relative error equal to in seconds, while it took seconds to obtain a relative error equal to with Algorithm 1. Note that for systems with a larger number of states the difference may be larger.

In Table 1, we compare the sample efficiency of Algorithms 1 and 2 with different input parameters, which we define as

where is the total number of generated samples, and is the number of samples in the set . Results in Table 1 indicate that Algorithm 2 has higher sample efficiency than Algorithm 1.

Our results also indicate that Algorithm 2 with , has on average a higher empirical convergence rate and a higher sample efficiency than Algorithm 2 with , . This indicates that combination of non-zero , improves convergence and sample efficiency, which can be explained as follows. When the total number of generated samples is low, we do not have sufficient information on the behaviour of the switching separatrix. Therefore we need to explore this behaviour by randomly generating samples, before we start minimising the relative error. This idea is similar to the so-called exploration/exploitation trade-off in reinforcement learning [28].

, , ,
, , ,
, , ,
, , ,
, , ,
, , ,
, , ,
, , ,
Table 2: Parameter values for systems in Subsection 5.2. The unspecified parameter values are the same as in (6).

5.2 Switching in the LacI-TetR System

Robust Switching in the LacI-TetR System. In Table 2, we specify different versions of the LacI-TetR system by varying some of the parameters in (6). After that we compute the switching separatrices and plot them in Figure 6. Note that the separatrices for , intersect, which happens since the vector fields and of systems