A simple electronic device to experiment with the Hopf bifurcation

TL;DR

This paper introduces a simple, low-cost electronic circuit that demonstrates the Hopf bifurcation, allowing for experimental study and control of dynamical behaviors in a physical system.

Contribution

The authors design and analyze a novel electronic circuit inspired by the genetic repressilator to experimentally demonstrate and control the Hopf bifurcation.

Findings

The circuit exhibits two distinct dynamical regimes with oscillations and steady states.

The bifurcation parameter can be easily tuned via a resistor value.

The system effectively models key features of the Hopf bifurcation.

Abstract

We present a simple low-cost electronic circuit that is able to show two different dynamical regimens with oscillations of voltages and with constant values of them. This device is designed as a negative feedback three-node network inspired in the genetic repressilator. The circuit's behavior is modeled by a system of differential equations which is studied in several different ways by applying the dynamical system formalism, making numerical simulations and constructing and measuring it experimentally. We find that the most important characteristics of the Hopf bifurcation can be found and controlled. Particularly, a resistor value plays the role of the bifurcation parameter, which can be easily varied experimentally. As a result, this system can be employed to introduce many aspects of a research in a real physical system and it enables us to study one of the most important kinds of bifurcation.

Tables

Table 1. Table 1: Frequencies as a function of R 3 subscript 𝑅 3 R_{3} for the experiemntal device.

$R_3$ [K$\Omega$]	Frequency [KHz]
1	1.88
2.2	1.38
5	1.14
10	1.11
20	1.08
47	1.07
100	0.872

Equations

$$\frac{d{V}_C}{dt}=\frac{V-V_C}{(R_1+R_2)C}.$$

$$\frac{d{V}_C}{dt}=-\frac{V_C}{R_{2}C}.$$

$$\frac{dV_{C}}{dt}=\frac{1}{\bigg{(}R_{1}u(V_{in})+R_{2}\bigg{)}C}\bigg{(}Vu(V_{in})-V_{C}\bigg{)}.$$

$$u(x)=\frac{1}{2}\Bigg{[}\tanh\Bigg{(}\beta\bigg{(}h-\frac{x}{R_{3}}\bigg{)}\Bigg{)}+1\Bigg{]}.$$

$$\beta =I_b/I_c$$

$$\frac{d{V}_i}{dt}$$

$$u(V_p)=\frac{V_p}{V}.$$

$$\epsilon(V_{p})=\Bigg{|}u(V_{p})-\frac{V_{p}}{V}\Bigg{|}.$$

$$A=\frac{\partial F_{ij}}{dV_{i}}\Bigg{|}_{\vec{P}}=-\frac{1}{\big{(}R_{1}u(V_{j})+R_{2}\big{)}C}\Bigg{|}_{\vec{P}}$$

$$\displaystyle B=\frac{\partial F_{ij}}{dV_{j}}\Bigg{|}_{\vec{P}}$$

$$u^{\prime}(x)=-\frac{1}{2}\frac{\beta}{R_{3}}sech^{2}\Bigg{(}\beta\bigg{(}h-\frac{x}{R_{3}}\bigg{)}\Bigg{)}.$$

$$J=\left[\begin{array}{ccc}A& 0& B\\ B& A& 0\\ 0& B& A\end{array}\right].$$

$${\lambda }_1$$

$${\lambda }_2$$

$${\lambda }_3$$

Full text

Hopf bifurcation, dynamical systems, electronic circuit, repressilator

I Introduction

Physical systems, generally, are mathematically modeled by differential equations and studied in the framework of the dynamical system theory Strogatz-book ; Perko . This formalism is also applied to many other systems like biological and economic ones Murray . The main goal in this theory is to find the behavior of the solutions or the trajectories in the space of system variables and, in particular, how they evolve to different kind of attractors, like fixed points and limit cycles. Notably, many systems have the property of having several qualitatively different solutions as a function of a certain parameter. The bifurcation analysis is the part of dynamical systems that studies the way in which these changes of dynamics arise with respect to these bifurcation parameters.

The Hopf bifurcation is one of the most important and known of the bifurcations in this theory. It is characterized by changes of dynamics between stationary states (stable fixed points) with fixed values for the system variables, and, oscillatory dynamics (limit cycles) where these variables evolve periodically in time. The system dynamics are controlled by a bifurcation parameter which is in general an important and characteristic parameter of the system.

This bifurcation emerges in many mathematical models that intent to describe different real systems. For example, the Hodgkin-Huxley HH-model and the FitzHugh-Nagumo FN-model models for the neural membrane potential. Non-linear chemical oscillators like the Belousov-Zhabotinsky reaction BZ-model with the Oreganator oreganator-model model as a possible mechanism of operation, and, the Brusselator brusselator-model . Some special cases of the predator-prey model PP-model for population dynamics. The Lorentz attractor lorentz-model as a paradigmatic example of deterministic chaos. And finally, we mention the repressilator repressilator-model as an example of a genetic network.

We can figure out from the previous paragraph that the Hopf bifurcation plays an important role in the characterization of many systems. However, most of the systems that the previous models intent to characterized are quite difficult to implement in a laboratory for undergraduate students. We propose in this work to design a very simple electronic circuit that is able to behave as having a Hopf bifurcation. We do not only observe typical oscillations, but we can also control the bifurcation parameter (resistor) in order to observe all the regimens of this bifurcation.

The proposed circuit is based on the genetic repressilator repressilator-model . It is an artificial genetic regulatory network that consists of three genes that repress each other in a loop structure. It is a paradigmatic example of the synthetic biology field synthetic-biology . Although very sophisticated electronic devices have been proposed for genetic systems paper-kurst , the main principle of operation of the repressilator can be mimicked with a simple $RC$ circuit. In effect, some of these circuits have been implemented in order to study synchronization properties on networks of artificial like-genes japoneses-circuitos .

This manuscript is organized as follows. In the next section, we introduce the electronic circuit and the mathematical model to describe it. We also develop a linear stability analysis and we find the eigenvalues of the stability matrix. In the third section, we present the analytical, numerical and the experimental results. Finally, we discuss this work and present the conclusion in the last section.

II Model

The genetic repressilator consists of three genes that repress each other in a loop-like structure as it is shown in Fig. 1(a). The flat arrows are the biological representation of a repressive directed interaction from a node to the other. We can give an insight of the operation of this system for certain parameter values as follows. Each node or gene is characterized to behave as a bistable unit with high and low levels of expressions. Thus, when a node has a high expression it represses strongly the node that it points to. The effect in the second node is to reduce its expression and takes the lower level. And so it results in a reduction of the repression of the third node that takes a high state. It is clear that this chain of effects cannot be stabilized since the system presents a kind of frustration and oscillatory dynamics emerge.

II.1 Circuit design

As a first task, we designed an electronic node that is able to present dynamics with two levels of expression. A node in this circuit can be seen in Fig. 1(b). This is a simple $RC$ circuit where the path of charge and discharge of the capacitor can vary depending on the transistor state which operates basically as a switch. The expression level of this device is given by $V_{C}$, and, the input signal is given by $V_{in}$.

In the first case, when $V_{in}=0$, the transistor is in the cut-off region and does not allow current to flow between its collector and emitter, and the circuit is described as

[TABLE]

Here, the capacitor will get a charge $V_{C}=V$ with a time constat $\tau_{1}=(R_{1}+R_{2})C$.

In the second case, when $V_{in}/R_{3}\gg h$, the transistor is in the saturation region and it conducts current between the collector and the emitter, and the system is described by

[TABLE]

In this case the capacitor discharges to zero with a time constant $\tau_{2}=R_{2}C$. The parameter $h$ indicates the minimum value of current needed in order to activate the transistor.

We note that when current $V_{in}/R_{3}$ in the base of the transistor is not enough to reach the saturation region the circuit behaves as a combination of the two previous cases. We can model this general situation like

[TABLE]

In this equation the step-like function

[TABLE]

The function $u(x)$ is an approximated description of the behavior of this circuit with a NPN transistor. Parameters $\beta$ and $h$ depend on the internal properties of the transistor and the other elements of the circuit. $\beta$ is the current gain of the transistor defined as:

[TABLE]

where $I_{b}$ is the base current and $I_{c}$ is the collector current. $h$ is a hybrid parameter of the transistor in common emitter configuration. It depends on the quiescent point of the transistor, i.e., the external components in the circuit. As a result of this design, the voltage $V_{C}$ can vary between $V$ and zero depending on $V_{in}$. In particular, $V_{in}$ plays the role of a repressive signal since $V_{C}$ behaves vice verse to $V_{in}$.

In the next step, we connect three of these nodes in a loop in order to construct a negative feedback closed-loop network. The output voltage $V_{i}$ of a node (previous $V_{C}$) is used as input ($V_{in}$) of the next node. The dynamics for such a system is as follows

[TABLE]

Here, only the elements $F_{12}(V_{1},V_{2})$, $F_{23}(V_{2},V_{3})$ and $F_{31}(V_{3},V_{1})$ are different from zero.

The model we obtain in (6) is a non-linear system of coupled first order differential equations. The non-linearity is the result of function $u(V)$ with its sigmoid shape ($\tanh()$). Additionally, as a result of this particular function, we cannot find a closed-form expression for the solution to this system. Note that this kind of expressions are common for modeling circuits with transistors. For example, the work cited in ref. japoneses-circuitos also employs a similar description for its device.

II.2 Stability analysis

We start this analysis of the mathematical model by searching the fixed points of the system (6), and later, we develop the linear stability analysis on this point.

II.2.1 Fixed point

The fixed point of the system (6) is the point in the space of voltages $\{V_{i}\}$ ($i=1,2,3$) where the condition $\frac{dV_{i}}{dt}=0$ is satisfied. In order to find this point we argue in the following way. The three coupled equations are identical among them with a strong symmetry with respect to the voltages (circular reference). As a result, they must fulfill the same condition simultaneously. That is:

[TABLE]

Here, $V_{p}$ is the voltage that the capacitors must have in order to settle the system on the the fixed point. Although we cannot find a close-form expression for $V_{p}$ we can observe that it corresponds to the intersection of a straight ($y=V_{p}/V$) and a negative hyperbolic tangent ($y=-tanh(V_{p})$). Thus, it is simple to see that there is only one possible intersection of these two curves. As a result the system has always only one fixed point. We call this point $\vec{P}=(V_{p},V_{p},V_{p})$.

Since we cannot find explicitly the value of $V_{p}$ because the characteristics of the function $u(x)$, we propose to get this value numerically by reducing the quantity

[TABLE]

This miminization is performed by variying $V_{p}$ between zero and $V$ each small intervals $\Delta V_{p}$ and computing the quantity $\epsilon(V_{p})$. A small value of $\Delta V_{p}$ ensures a good aproximation for the position of the fixed point.

II.2.2 Linear stability analysis

We evaluate the Jacobian matrix of the system (6) on the fixed point $\vec{P}$. The elements of this matrix are

[TABLE]

and

[TABLE]

with

[TABLE]

The Jacobian matrix has the shape

[TABLE]

This matrix has the following eigenvalues

[TABLE]

In order to know the kind of dynamics as a function of $R_{3}$, we need to find the sign of the real parts of these eigenvalues. In particular, we need to know if $Re(\lambda_{2,3})=0$ for a critical value $R_{c}$ for the bifurcation parameter $R_{3}$. Unfortunately, we cannot find an analytical formula for this dependence and we need to evaluate these expressions numerically.

III Results

In this section, we present several results of this device and its mathematical model. We show numerical simulations and experimental work.

III.1 Experimental study

We constructed the experimental device whose design is shown in Fig. 1(b). The values of the electronic elements are: $R_{1}=1K\Omega$, $R_{2}=1K\Omega$, $C=220nF$, $V=5.45V$. The transistor is type $NPN$ model $2sc2235$. The experimental setup is shown in Fig. 2. We note that in the experimental device there are variations on the values of the electronic components between nodes because of the manufacturing dispersion. As a result, the analytical derivations we have performed before are qualitatively correct if these variations are relatively small.

We measure the voltages $V_{i}$ in the capacitors using an oscilloscope as a function of different values of $R_{3}$. The initial condition of the device is discharged with the power supply off. After turning on the voltage source $V=5.45V$ and after a transient the output voltages $V_{i}$ oscillate or reach constant values in time. Fig. 3 presents four measurements for different values of $R_{3}$ where oscillations are established. We find that the frequencies are higher for smaller values of $R_{3}$ and that oscillations cannot be found after a critical value $R_{c}=110K\Omega$ for $R_{3}$ (Fig. 3(d)). Table 1 summarizes the set of measurements for this device.

We observe in Fig. 3(a) that the voltages of the capacitors possess an exponential charge and discharge that can be easily seen. In effect, with a frequency of $1.88KHz$ the period $T=5.32\times 10^{-4}s$ is around the addition of the two characteristic times of this system: $\tau_{1}=R_{2}C=2.2\times 10^{-4}s$ and $\tau_{2}=(R_{1}+R_{2})C=4.4\times 10^{-4}s$. On the other hand, when the system is close to the bifurcation point (Fig. 3(c) with $R_{3}=100K\Omega$) the signals are almost sinusoidal as we expect from the theory of Hopf bifurcations.

Finally, we have measure the three output voltages when oscillations are absent for $R_{3}\gg R_{c}$ ($R_{3}=200K\Omega$). We find that the voltages are: $V_{1}=2.5V$, $V_{2}=2.93V$ and $V_{3}=2.96V$. Note that the three voltages are different, and this differs from the theoretical analysis from which we expect the same values. This variations in the fixed point position are in part due to the dispersion of values of the electronic components.

III.2 Numerical analysis

In this subsection we present numerical studies in order to compute the voltage values on the fixed point and the eigenvalues of the stability matrix. We also share some examples of the possible dynamics that this model can show. In order to have compatible results in our simulations with the experiment, we chose $\beta$ and $h$ values that make fit both behaviors in the region where the Hopf bifurcation take place. As a result we take $\beta=8.5\times 10^{4}$ and $h=1.8\times 10^{-5}$.

Fig. 4(a) presents the position of the fixed point $V_{p}$ as a function of the bifurcation parameter $R_{3}$. The solid curve is calculated numerically by reducing the quantity $\epsilon(V_{p})$ (eq. 8). We also include the experimental values obtained for $R_{3}=200K\Omega$ (blue dots). We observe that these points are similar qualitatively to the proposed model.

Fig. 4(b) and 4(c) present the real and imaginary parts of the three eigenvalues of the Jacobian matrix as a function of $R_{3}$. We find that $\lambda_{1}$ is real and always negative indicating that the fixed point is stable in its direction. The eigenvalues $\lambda_{2}$ and $\lambda_{3}$ are complex conjugates. Their $Re(\lambda_{2,3})$ is positive for values smaller than a critical value $R_{c}=112K\Omega$ and negative for larger ones. This indicates that the fixed point is unstable for $R_{3}<R_{c}$ and stable for $R_{3}>R_{c}$. In $R_{3}=R_{c}$ we have the Hopf bifurcation point. In this point $|Im(\lambda_{2,3}(R_{c}))|=5582$, thus, the frequency of the oscillations close to the onset of the bifurcation is $888Hz$. Note that with our election of $\beta$ and $h$ we can fit the main properties of the Hopf bifurcations with our model.

We present some characteristic trajectories in Fig. 5. They are performed by integrating the system (6) with a Runge-Kutta algorithm of second order and a time step $\Delta t=1\times 10^{-5}s$. In Fig. 5(a) we show an evolution far from the onset of the bifurcation ($R_{3}=1K\Omega$). We observe that the voltages show that the capacitors charge and discharge almost completely, generating signals with big amplitudes. The frequency of this signals is $317Hz$. Note that these results are quite different from our experiment where the device has a frequency of $1.88KHz$ and smaller amplitude in the signals. These differences arise because our model does not represent well the transistor in the linear region of operation. The transition between cut off and saturation region is determined by the resistor $R_{1}$. In this circuit the transistor enters to the active region before the capacitor completely discharges.

Fig. 5(b) presents the voltages as a function of time close to the bifurcation point with $R_{3}=100K\Omega$. We observe that the signals have a sinusoidal shape and relatively small amplitude. These two characteristics are typical of the Hopf bifurcation. The frequency of these signals is approximately $910Hz$ which are in the order of the expected value. In Fig. 5(c) we show the evolution for $R_{3}=200K\Omega$ where the fixed point is stable. We observe that all the signals decay with oscillations to a constant value.

Finally, in Fig. 6 we present the three previous trajectories in the space of voltages $\{V_{i}\}$. All the trajectories start from the same point and they evolve in quite different ways. When $R_{3}=200K\Omega$ the trajectory is a spiral going to the fixed point (black curve). Close to the bifurcation point ($R_{3}=100K\Omega$) the trajectory evolves to a limit cycle (red curve) with an almost circle shape on a plane. In the last case, with $R_{3}=1K\Omega$, the limit cycle has a shape far from a circle and it is no longer contained in a plane (blue curve).

IV Discussion and conclusions

In this work we have presented a simple and low-cost device with a rich behavior that can be modeled by a dynamical system with a Hopf bifurcation. This circuit does not require a deep understanding of electronic and can be handled by undergraduate students that have taken the basic physics courses. Notably this circuit can be seen as a physical implementation of the negative feedback three-node motif that is intensively study in biology Milo and conceptually similar to the genetic repressilator.

The performed analyses allow to understand different aspects of the proposed model. In effect, with the linear stability study we can find the the critical value of $R_{3}$ and estimate the frequency of the system close to the bifurcation point in an almost analytical way. From this analysis we also observe that the resistor $R_{3}$ plays the main role as bifurcation parameter. This is because it enters as argument of the function $u(V)$ that is the kernel of the non-linearity of this system and it is the responsible to switch between the two possible states of a node. On the other hand, resistors $R_{1}$ and $R_{2}$ have the role of controlling the characteristic times for charge and discharge of the capacitors and they are not a source of non-linearity for the system.

Numerical integration of the model allows us to find the trajectories of the system in the space of voltages. We can use this information to compute the main characteristic of the signal like frequency and its shape for the whole set of possible values of the system parameters. Note that these values cannot be analytically determined since the system (6) does not allow a close-form expression for its solutions.

In conclusion, this system allows us to study important concepts of dynamical systems such as stability analysis and bifurcation theory. It allows also to perform analytical and numerical studies of the proposed model. Finally, we can construct a device where the theoretical results can be contrasted with the measurements.

PK acknowledges financial support from SeCTyP-UNCuyo (project M028 2017-2018) and from CONICET (PIP 11220150100013), Argentina. DNR acknowledges financial support from EVC-CIN fellowship 2016, Argentina.