Title image Project #3

Project 3: Simulation of coupled circadian oscillators (due Friday Mar 2 at midnight)

This week, we will incorporate intercellular signaling (also called "coupling") into our model of the circadian clock. If all goes well, then we will see that oscillators with different intrinsic periods can be cajoled into oscillating together. We will simulate it with increasing coupling strengths to determine which coupling strength works best for a particular distribution of periods. We will also measure the performance by timing our runs.

  1. One Simulation at a Time, with Coupling

    The first task is to incorporate support for VIP signaling. This means, we will need to add a parameter and update the code that evaluates the phase velocity. At every time step, we will need to determine the average VIP output and use it to determine how much the VRC affects the phase velocity. Be sure to read Stephanie's Reference Guide before you begin. This time, focus on the VRC and the mean field coupling equations.

    1. Update the definition of the Params struct in phase_support.h to include a parameter that will control the strength of the VIP output signal:
      typedef struct{
        float *periods;
        float vipStrength;
      } Params;
    2. Update sim_slow.c to include a VIP strength parameter.
        if (argc < 2) {
          printf("Usage: ./sim_slow  [ []] \n");
          printf("  where\n");
          printf("        id's the file which stores the output\n");
          printf("            Use the .phs extension in this filename\n");
          printf("        is the number of oscillators to simulate\n");
          printf("         (defaults to 100)\n");
          printf("        is the floating point strength of VIP\n");
          printf("            (default value is 0.0)\n");
          return 1;
        char *filename = argv[1];
        int Nx = 100;
        if (argc > 2) {
          Nx = atoi( argv[2] );
        Params params;
        if (argc > 3) {
          params.vipStrength = atof( argv[3] );
          printf( "setting vip strength %f\n", params.vipStrength );
        else {
          params.vipStrength = 0.0;
    3. Update phase_support.c to implement VIP output and to include the VRC. It is up to you how to design the code.
    4. Test it by running sim_slow. As in project 2, you can compare Stephanie's .phs files with ones you create. The first call should reproduce the data from last week. The second two add coupling. Here are the commands used to create the files.
      ./sim_slow uncoupled_5.phs 5 0.0
      ./sim_slow coupled_5_1.phs 5 0.1
      ./sim_slow coupled_5_10.phs 5 1.0

      You can access the new files from the web:

      • curl -O http://cs.colby.edu/courses/S18/cs336/projects/proj03/proj03/coupled_5_1.phs
      • curl -O http://cs.colby.edu/courses/S18/cs336/projects/proj03/proj03/coupled_5_10.phs
  2. Lots of Simulations for one VIP strength

    Ultimately, our goal is to run many simulations for each VIP strength starting from random initial conditions. First, let's write code that runs multiple simulations from random intitial conditions for one VIP strength.

    1. Add support for random initial conditions.
      1. Add code to phase_support.h:
        // Return a float array with Nx periods taken from a Gaussian distribution
        // with the given mean and standard deviation.
        float *generateGaussianPeriods( int Nx, float mean, float std_dev );
      2. Add code to phase_support.c:
        // Adapted from the C++ code on Wikipedia
        // https://en.wikipedia.org/wiki/Box–Muller_transform
        double generateGaussian(double mu, double sigma)
        	static const double epsilon = DBL_MIN;
        	static const double two_pi = 2.0*3.14159265358979323846;
        	double z1, u1, u2;
        	   u1 = rand() * (1.0 / RAND_MAX);
        	   u2 = rand() * (1.0 / RAND_MAX);
        	while ( u1 <= epsilon );
        	double z0;
        	z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2);
        	z1 = sqrt(-2.0 * log(u1)) * sin(two_pi * u2);
        	return z0 * sigma + mu;
        // Return a float array with Nx periods taken from a Gaussian distribution
        // with the given mean and standard deviation.
        float *generateGaussianPeriods( int Nx, float mean, float std_dev ) {
            // Students write this code. 
      3. Test generateGaussianPeriods by writing a new main program or by altering one of the existing ones.
      4. Download the main program: curl -O http://cs.colby.edu/courses/S18/cs336/projects/proj03/proj03/sim_stats_1s.c
      5. Run stim_stats_1s. Here are example output from Stephanie. Yours do not need to be identical, but should have answers in the same ballpark:
        ./sim_stats_1s 10 100 0.5 5
         1.112769 1.135372 1.090012 1.164543 1.060384
        Ran in 0.156025 seconds
        ./sim_stats_1s 10 100 0 5
         4.767481 4.829017 4.753807 5.007208 4.616753
        Ran in 0.170441 seconds
    2. Lots of Simulations for multiple VIP strengths

      1. Download the main program: curl -O http://cs.colby.edu/courses/S18/cs336/projects/proj03/proj03/sim_stats_ns.c
      2. Run sim_stats_ns to verify it is working. Shown below is the output when Stephanie runs it for 10 cycles with 100 oscillators and 5 trials per VIP strength:
        ./sim_stats_ns 10 100 5
        VIP strength 0.000000:  4.767481 4.829017 4.753807 5.007208 4.616753
        VIP strength 0.100000:  3.480247 3.362637 3.356769 3.363106 3.149650
        VIP strength 0.200000:  2.428410 2.397236 2.487279 2.523321 2.568427
        VIP strength 0.300000:  2.176871 2.010643 1.728647 1.746991 1.769865
        VIP strength 0.400000:  1.714436 1.481059 1.349530 1.314250 1.316329
        VIP strength 0.500000:  1.087567 1.076648 1.053374 0.945921 1.083624
        VIP strength 0.600000:  0.975974 0.833561 0.914562 0.850786 0.940000
        VIP strength 0.700000:  0.803137 0.823165 0.655570 0.788306 0.726310
        VIP strength 0.800000:  0.703527 0.654615 0.615664 0.628779 0.634605
        VIP strength 0.900000:  0.541345 0.626822 0.533902 0.604878 0.602841
        VIP strength 1.000000:  0.491972 0.540740 0.490722 0.477447 0.585429
        VIP strength 1.100000:  0.423047 0.476743 0.424007 0.442713 0.503413
        VIP strength 1.200000:  0.426519 0.471773 0.408094 0.455949 0.450449
        VIP strength 1.300000:  0.391024 0.391191 0.354587 0.343262 0.382828
        VIP strength 1.400000:  0.333713 0.429654 0.335122 0.337064 0.346848
        VIP strength 1.500000:  0.305867 0.294562 0.320573 0.301878 0.332879
        VIP strength 1.600000:  0.318930 0.303940 0.337960 0.298623 0.264145
        VIP strength 1.700000:  0.284260 0.259595 0.254676 0.320787 0.310911
        VIP strength 1.800000:  0.283433 0.266384 0.273583 0.276506 0.246806
        VIP strength 1.900000:  0.251058 0.273244 0.284011 0.253571 0.266674
        VIP strength 2.000000:  0.000000 0.000000 0.000000 0.000000 0.000000
        Ran in 3.580201 seconds
        Note that at a VIP strength of 2.0, it looks like everything is perfectly in sync. And, in a way, it is. But if you run the simulation with VIP strength of 2.0 using sim_slow, you will find that it halts all the oscillators.
      3. Run sim_stats_ns for 3 different period distribution widths (e.g. choose a standard deviation of 0.1, 0.25 and 0.5 h). Determine which VIP strength works best (allows for the most similar set of CT6-crossings in the 10th cycle) for each one. Do we need more coupling to make up for a wider range of periods? You should run sim_events and/or sim_slow to verify that the simulations are actually oscillating.


    • The relationship between the timing of the signal and the VRC is critical. Right now, the VIP signal is sent around CT6, which is also when the VRC crosses zero with a negative slope. It is consistent with the theory of entrainment that such a relationship should allow the oscillators to synchronize (mutually entrain). An interesting set of simulations would adjust the timing of the VIP signal to see how much worse it gets as it moves away from this ideal relationship. Another interesting idea for your write-up is to read up on the theory of entrainment and explain why this relationship works.
    • Expand your study in the final task to consider a wider variety of distributions.
    • Use an explicit network topology instead of a mean field (see Stephanie's Reference Guide for more information).
    • Expand your study of the timing performance of sim_stats_ns.

    Writeup and Handin

    1. Create a file named README.txt. Describe the code files (the purpose of each file) along with how to compile and run the code.
    2. Create a second file for your project report. This one should allow for images or tables, so should be in a format such as .docx or .pdf. This week, I want you to talk about (a) what was the most challenging aspect of the project, (b) which extensions you did, (c) how long it takes to run sim_stats_ns for a range of values of N (the larger the N, the better), and (d) the results from your final task that attempts to determine the best VIP strength for 3 period distributions. Spend 1-2 hours on this part (writing and running together).
    3. You should hand in all code necessary to run your solutions. Place all necessary .h, .c, and Makefile files in proj03 directory along with the README file. Stephanie will probably want to compile and run the code. It should be possible to do so without looking for any more files. Tar/zip up the directory and put it in your Private folder on Courses/CS336.