Quick start guide for molecular circadian clock model
GonzeGoldbeter model[1] †For example, let us consider the circadian oscillator model, which has been proposed by Goldbeter and the coworkers, with an external periodic force in Neurospora. The external periodic force corresponds to a lightdark cycle. Then, the circadian oscillator can be entrained to the external periodic force under a certain condition. However, when the period of the external force and the amplitude of one are changed, the NeimarkSacker bifurcation occurs. Let us calculate the NeimarkSacker bifurcation set in the plane. The circadian oscillator model is described by: In the following, we make the following programs, and use these:
Configurations †
Creation of the Programs †The following steps are done in a subdirectory. Any name can be used for the subdirectory. First, the following is executed in the MATLAB command window: >> BUNKI After the execution, a window named Select BUNKI Project Directory will open. Here the Project to be analyzed is created or selected. At this stage, no Project is generated. Therefore, input a Project name for identification in the Name field and click the OK button. In this case, since the circadian oscillator model is proposed by Goldbeter and is analyzed in detailed by Gonze, GG is used as the project name.*1 When creating a new Project, you are asked whether it is okay to create a directory (or folder). Please select Yes. SE (system editor) will start. We now create the programs for analysis.*2
PP †Here, we import the initial values of parameters and state variables to draw a phase portrait and to capture a stable periodic solution. First, start the PP (phase portrait) program from SE. Select Program > PP. The drawing window and the control panel will start. In the first run, the drawing window is the (x[1],x[2]) plane.*4 Here, there is no need to change the drawing window. Next, we deal with the control panel. First, we input the values of Initial states and Parameters as the initial condition.*5. The initial points to start the simulations are input (x[1],x[2]) in the Initial states field. Here, let us start from the origin (0,0) ! Please do not forget to input the parameter values before drawing. When first ran, 1 is entered for all parameters. Parameters are configured as follows.
Next, click the Start button to start simulations. Can you see that the orbit of a solution starting from the origin converges to a periodic solution ? The blue line and red point indicate the orbit of the periodic solution and the points of the stroboscopic map every time of , i.e., the point of the Poincare map. Parameters can be changed in real time. Let us change the parameter omega. The parameter Omega is related to the period of the sinusoidal force. Currently, the period of this external force is set to about 24 [h](T = 2π/omega) . Since the default step size to change the parameter is 0.1, which is too large, first we make changes here. Select menu: Setting > Parameter. Then, the configuration window will open.*6 Change the step size for omega displayed at the bottom from 0.1 to 0.001. It is OK after the change is completed. Now, what happens when the parameter omega is increased ? By increasing the parameter omega, the stable periodic solution changes to the unstable one. For instance, the following figure shows an attractor observed in the system when omega = 0.317. We make this figure more easier to understand. Uncheck the check box Orbit in Plot field in main panel. The orbit of the solution will not be displayed, and only a point sequence is displayed in the phase portrait. As time advances, the sequence of points settle on a closed curve. This is called the invariant closed curve (ICC). The emergence of the ICC is involved in the generation of the NeimarkSacker bifurcation.*7 In this case, by increasing the value of omega, the stable periodic solution changes to unstable one, and then the NeimarkSacker bifurcation (supercritical NS bifurcation, which a stable ICC can be observed.) occurs.~
Here, we return to the original settings by reducing omega. Then, it can be achived to click the button to decrease the parameter or to directly change the value of the parameter to its origibal value. You can see that the orbit of the solution converges to a stable periodic solution as expected. Next, we check the check box states flow in the Plot field. In this way, we can check the current values of the state variables in the list box during numerical integration. The values of the state variables sampled once per 2π/omegaperiod is output in the list box. In this example, we can see to converge to a stable fixed point of the Poincare map since the same number repeatedly appears in the list box. Here, by clicking the Clear button, the orbit during the transient response shown in the drawing window can be erased, and the steady state can be observed. Uncheck the check box Orbit, and stop drawing the orbit of the solution. Then only one point corresponding to the stable fixed point is shown on the drawing window since the system converges to a stable fixed point. Next, let us investigate bifurcations of this periodic solution (the fixed point of the Poincare map). The eigenvalues of the variational equation for the fixed point are evaluated when the parameters are changed. For this purpose, the FIX tool can be used. We pass the initial values to use FIX from PP. The important thing is that the values of the state variables and parameters at the steady state must be passed to FIX. To confirm whether the system reached a steady state, either check the check box states flow switch and confirm that the values of the state variables no longer change as described above, or click the Status button a few times and confirm that the information displayed in the list box does not change. After confirming that there are no more changes, select Export current status from Tools in the menu bar. With this procedure, the values of state variables, the parameter values, and other configuration values at the steady state are preserved and can be passed to the next program (see section on FIX). FIX †Here, the objective is to track the accurate location of a periodic solution when the parameters are changed and also to monitor the eigenvalues during tracking. Users must adjust the configurations to run FIX. For simplicity, here let us run FIX as is. Normally, the FIX program is run after the values saved by the PP program (using Export Current status) is imported as the initial values for Newton’s method.*8 First, let us start up the FIX program from the PP control panel. Select Program > FIX. When initially starting up the FIX, default values are entered in each field. Information on the fixed point for a periodic solution obtained with PP is imported as the initial configuration value. Select Tools >Import initial point. After the selection, the information obtained with PP such as the parameter values and the coordinates of the fixed point for the periodic solution are reflected. By using PP, we saw that the dynamics of the system changes by increasing the parameter omega. Therefore, select the parameter to be changed as omega and configure the step size to increase the parameter. Click the Start button to start calculations. A new window will open when the calculations start. This window is to monitor the location of the eigenvalues on a Gauss plane. The NeimarkSacker bifurcation of a periodic solution occurs when the absolute value of the complex conjugate eigenvalues becomes "1", that is when passing through the unit circle. We can observe that as the parameter changes, the location of the eigenvalues change, and the pair of eigenvalues on the window approaches "unit circle". Parameters are continuously changed, and when a bifurcation point is detected, the program stops.*9 In this example, the NeimarkSacker bifurcation is detected because the absolute value of the complex conjugate eigenvalues became 1. The values of the parameters and the state variables when the program stopped can be used as the initial point for BF, the program to track bifurcation curves. The information when the bifurcation point is detected are exported for BF. When Export BF point is selected here, the status, parameter values, and other necessary configuration values are preserved. The preserved values are passed to the next program (see section on BF). BF †The BF tool is used to obtain a bifurcation curve by changing system parameters. The user must adjust the configuration to run BF. For simplicity, here let us run BF as is. Normally, the BF program is run after the values of the state variables and parameters at the bifurcation point detected by the FIX program are imported as the initial values for Newton's method.*10 First, let us run the BF program from the FIX main panel. Select Program > BF. When starting up BF at the first time, default values are entered in each field. Various information on the bifurcation point obtained with FIX is imported as the initial value. Select Tools >Import BF point from main panel. After the selection, various information such as the type of bifurcation, the period of the periodic solution, and the angle of eigenvalues are input automatically into each field. To tracke the bifurcation set on the (omega,vsmax) plane, check the radio button omega as the X parameter and vsmax as the Y parameter. Furthermore, configure items such as step size to change the parameters and the value to stop the BF when reaching end parameter. To save the calculated results as a data file, input the output file name in Output. After these configurations, click Start button to start tracking of the bifurcation curve. Then, the Eigenvalues window and Bifurcation diagram window will open as new windows. The Eigenvalues window is used to visualize the change in the root locus of the eigenvalues during calculation of the bifurcation curve. Furthermore, to monitor the change of the root locus, we can visually understand occurrence of a bifurcation point with codimension 2. The Bifurcation diagram window is a window to monitor how the NeimarkSacker bifurcation is being tracked in real time. (b)：Bifurcation diagram in the (omega,vsmax)plane. The calculation result can be monitored in real time.
