Design Principles of the Yeast G1/S Switch

Single-cell microscopy and computational modeling offer novel mechanistic insight into the G1/S switch that initiates DNA replication in budding yeast, revealing a Clb5/6-Cdk1 and Sic1 feedback loop and new rules of biochemical circuit design.

Here, lowercase symbols, kinase1, kinase2, and inhibitor denote genes; whereas uppercase symbols, KINASE1, KINASE2, and INHIBITORp i denote proteins. INH_K2p i represents the Inhibitor--Kinase2 complex. Rate of each reaction is given after the comma. The reactions for KINASE1 and KINASE2 synthesis employ Hill functions to produce sigmoidal curves. These reactions are implemented as a function of time. The reaction rate is calculated using the function Hill , , , = represents time in seconds. The indices, i and j are integers denoting the number of phosphorylated sites (0 ≤ i ≤ 5 and 0 ≤ j ≤ 6). In the 1P model, i=0 and j=0, 1. We used molecule counts in the simulations (not concentrations), and chose k s1 =k s2 =160, k d1 =k d2 =0.08, K s1 =K s2 =10, n s1 =n s2 =4, k s3 =23.1, k 3 =0.5, k 4 =0.5, k 5 =0.24, k 6,i<6 =0.0154, k 6,6 =0.17. Constants, k s1 , k s2 , k d1 , k d2 , K s1 , K s2 , n s1 , and n s2 are subject to extrinsic noise (to simulate cell to cell variability). Each time a copy of the circuit is run, the specified reaction rate is multiplied by 1+ρ, where ρ is a random number uniformly distributed on (-x, x). x is chosen as 0.2 for k s1 , k s2 , K s1 , and K s2 ; and 0.7 for n s1 and n s2 . The phosphorylation rates, k 1,i , and k 2,i are real numbers between [10 -10 , 0.1] and are subject to mutations. At the beginning of the simulation, they are assigned random numbers from a log-uniform distribution on [10 -6 , 0.1]. Dephosphorylation was assumed to occur at a constant rate. Mutations are modeled by multiplying the reaction rates by a random number uniformly distributed on (0, 2]. Mutation rate per reaction rate per generation was chosen as 0.3. Mutations are carried from one generation to the next, unlike the extrinsic noise, which is not heritable. k 1,i = 0 in the double-negative feedback loop simulations, and k 2,i = 0 in the linear circuit.
We define timing of the Kinase2 activation as the time the number of Kinase2 exceeds the number of Inhibitor. We take sharpness as the slope of the free Kinase2 curve at its half-maximum level.

Optimization by in silico Evolution
The algorithm resembles evolution by natural selection, and works as follows.
1. Initialization. N=1000 copies of the network is generated. Phosphorylation rates of each copy (realization) are assigned random rates as discussed above. Initial protein counts are all zero except Inhibitor=1500. All genes have exactly one copy each.
2. Scoring. Reaction rates that are subject to extrinsic noise are modified as described above. (In this set, only Kinase1 and Kinase2 production/degradation are subject to noise.) Then, the network is run, and a fitness score is calculated based on the output. Fitness score for timing is Here t desired denotes desired activation time for Kinase2. Sharpness score is simply the sharpness of the free Kinase2 curve as defined above. Lower scores mean higher fitness.

Elimination and duplication.
Half of the population with lower fitness is eliminated. The rest is duplicated to make up the deficit. 4. Mutation. Duplicates are mutated as described above. Mutations only affect phosphorylation rates (i.e. catalytic efficiencies).
Steps 2-4 are repeated for m=1000 generations. Longer runs were not observed to change the final distributions of the phosphorylation rates qualitatively.