SUPPLEMENTARY NOTEBOOK
Elucidating the Genotype-Phenotype Map by Automatic Enumeration and Analysis of the Phenotypic Repertoire
Jason G. Lomnitz\(^1\) and Michael A. Savageau\(^{1,2,*}\)
\(^1\)Department of Biomedical Engineering, and \(^2\)Microbiology Graduate Group, University of California, Davis, CA 95616 USA
\(^*\)To whom correspondence should be addressed. Tel: +01 530 754 7350; Fax: +01 530 754 5739; Email: masavageau@ucdavis.edu
Classification: BIOLOGICAL SCIENCES – Biophysics and Computation Biology
Keywords: System design space, circuit architecture, mode of transcription control, genetic oscillators, systems and synthetic biology, model deconstruction and reduction
Please note: This analysis uses v0.2.1 of the Design Space Toolbox V2 and Python interface. These versions can be downloaded using the toolbox update script.
./toolbox_update.py -t 0.2.1 -i 0.2.1IPython Notebook Showing Automated Analysis of Chemical and Biochemical Systems
The following IPython notebook reproduces the results published in the article entitled "Elucidating the Genotype-Phenotype Map by Automatic Enumeration and Analysis of the Phenotypic Repertoire". This notebook is structured to mirror the Online Methods section in the Supplementary Materials.
First, we analyze a small circuit involving two regulators creating a positive feedback loop (Eqs. S1-S2 in the Online Methods). Then, we apply the same methods to more complicated circuitry that includes a relaxation oscillator (Eqs. S35-S38 in Online Methods) and analyze a general class of two-gene circuits with several genetic oscillator designs (Eqs. S45-S50 in Online Methods).
Design Space Toolbox V2
The analysis shown here uses the Design Space Toolbox V2 (work in progress) for Python. The Design Space Toolbox V2 is a C library (libdesignspace) for the analysis of biochemical systems through a method of phenotypic deconstruction (1-3). This tool is based on the Design Space Toolbox 1.0 for Matlab® (4), with many additions that improve the automated construction and the automated analysis of biochemical systems. The library is accesed within an interactive python environment that facillitates rapid cycles between model development, model analysis and critical hypothesis testing.
Analysis with the Design Space Toolbox V2 is done through the dspace package that must be imported into the python environment. This package automatically adds the most commonly used classes for model analysis to the dspace namespace,
import dspace
The dspace.plotutils module is imported to add methods for data visualization,
import dspace.plotutils
by importing the dspace.plotutils module we add plotting methods to classes in the dspace package. The plots are created through matplotlib (5), with dependencies on numpy and scipy (6) using the IPython notebook environment (7). This combination of packages provides a robust platform for scientific computing (8,9).
import matplotlib
import numpy
import scipy
# We create an alias for the matplotlib.pyplot module
import matplotlib.pyplot as pyplot
# We create an alias for the numpy package
import numpy as np
We import the dspace.display module to add rich display of Equations within the notebook environment.
import dspace.display
For reproducibility, the work shown in this notebook has been generated with the following versions of these packages,
print 'libdesignspace ' + dspace.__c_toolbox__.__version__
print 'dspace ' + dspace.__version__
print 'matplotlib ' + matplotlib.__version__
print 'numpy ' + numpy.__version__
print 'sciPy ' + scipy.__version__
Once these packages have been imported into the python environment, we proceed with the analysis of mathematical models of biochemical systems using the dspace package. Model analysis is done through the following series of steps: we (a) define a mathematical model using a string representation, (b) instantiate the Equations object and (c) instatiate the DesignSpace object. Once the DesignSpace object is created, analysis of the model can be directed through a variety of methods to test specific hypotheses. In the next section, we use the simple system defined by Eqs. S3-S6 in Online Methods to illustrate the methods that will be used throughout this notebook. Then, we apply these methods to more complex circuitry that can produce complex behaviors such as hysteresis and oscillations.
Automated Enumeration of Qualitatively Distinct Phenotypes
The aim of this section is two-fold: first, we show a general overview of some of the available methods for analysis of models of complex biochemical systems; and second, we familiarize the reader with the methods that will be used in the analysis of a class of gene regulatory circuits. The system used throughout this section is defined by the recast Eqs. S3-S6 in Online Methods. This system has a relatively simple design involving two dependent variables, one independent variable and 8 parameters. However, despite its size it is capable of displaying a variety of phenotypes with interesting behaviors.
For a model to be analyzed using the dspace package it must be a system of ordinary differential equations and algebraic constraints in the generalized mass action form. Each equations is represented by a string that is parsed to create an internal representation of the system. Differential equations are defined using Newton notation such that X1. \(= dX_1/dt\). For example, Eq. S3 in Online Methods is represented by the following string,
'X1. = a1*X4^-1 + a1*rho1*X2^2*K1^-2*X4^-1 - b1*X1'.
Similarly, algebraic constraints are defined with the left-hand side set equal to \(0\), and the corresponding auxiliary variable is explicitly indicated during creation of the Equations object. Alternatively, the left-hand side of the algebraic constraint for the auxiliary variable can be explicitly defined.
The model for this section is represented by the following list of strings,
equations = ['X1. = a1*X4^-1 + a1*rho1*X2^2*K1^-2*X4^-1 - b1*X1',
'X2. = a2*X1*X5^-1 + a2*X1*rho2*X3^2*K2^-2*X5^-1 - b2*X2',
'X4 = 1 + X2^2*K1^-2',
'X5 = 1 + X3^2*K2^-2'
]
The string representation of the model uses a simplified notation for the parameters. For example, \(\alpha_1\) in the model is represented by the string a1. However, the simplified representation can be converted automatically by the dspace package using simple symbol substitution. Therefore, to improve the display of the equations and parameters, we specify a dictionary that relates the simplified notation to symbols using the latex format.
symbols = {'X1':'X_1',
'X2':'X_2',
'X3':'X_3',
'X4':'X_4',
'X5':'X_5',
'a1':r'\alpha_1',
'a2':r'\alpha_2',
'b1':r'\beta_1',
'b2':r'\beta_2',
'rho1':r'\rho_1',
'rho2':r'\rho_2',
'K1':'K_1',
'K2':'K_2'}
The model in string representation is the primary input into the dspace package. We use the list of equations input in In [6] to create an Equations object.
equations_object = dspace.Equations(equations, latex_symbols=symbols)
The equations for this model using the symbols we have defined are displayed below,
equations_object
We use the Equations to instantiate a DesignSpace object that is responsible for automatically deconstructing a complex model into qualitatively distinct phenotypes, and reassembling the results of the deconstruction.
ds = dspace.DesignSpace(equations_object)
A very important component of the phenotypic deconstruction is identifying the upper bound on the number of potential qualitatively distinct phenotypes for a particular system. This is calculated based on the number of positive and negative terms in each of the system equations. The number of positive and negative terms in each equations can be shown as a sequence of numbers, that is collectively referred to as the design space term signature. The term signature can be readily obtained from the DesignSpace object.
print 'Term signature: ' + ds.signature
The maximum number of potential phenotypes is given by the number of combinations of dominant positive and dominant negative terms for each equation, and it is equal to the product of all of the elements in the term signature (1). This is calculated internally and the information can be retrieved from the DesignSpace object,
ds.number_of_cases
Thus far we have shown the basic steps for creating the DesignSpace object. In general, we proceed with the analysis of the model from a global perspective that then directs more focused analyses to specific regions of parameter space (10). However, here we show the analysis of individual qualitatively distinct phenotypes first, to make the concept of qualitatively distinct phenotypes concrete and to reveal the implications of this phenotypic deconstruction. Detailed definitions of these concepts are found in the Online Methods section of the article and elsewhere (1, 4).
Briefly, a unique combination of dominant positive and dominant negative terms for each of the equations of the system defines a potentially valid, qualitatively distinct phenoype of the system that is uniquely identified by its case number or its case signature. For example, the phenotype identified as Case 2, with Signature [11111121], is obtained from the DesignSpace object in the form of a Case object. The Case object can be inspected for its case number, case signature, conditions for term dominance, amongst many other properties.
case2 = ds('11111121', by_signature=True, constraints=['rho1 > 1', 'rho2 > 1'])
print 'Case Number: ' + str(case2.case_number)
print 'Case Signature: ' + '[' + case2.signature + ']'
print 'Dominance Conditions:'
case2.conditions
A qualitatively distinct phenotype is represented by a unique sub-system (S-System). The S-System corresponding to each Case object is contaned in an SSystem object that is readily available. In the following block we obtain the SSystem for Case 2 and display the S-System equations.
ssystem = case2.ssystem
print 'S-System Equations:'
ssystem.equations
The steady-state solution for S-System models are calculated automatically by the dspace package, and these can be displayed,
print 'Solution:'
ssystem.solution
The steady-state solution of S-System equations can be calculated analytically by virtue of the fact that the steady-state solution is linear in logarithmic space. In the block below, we show the steady-state solution for the S-System of Case 2 in logarithmic coordinates.
print 'Solution (logarithmic):'
ssystem.solution_log
If an S-System has a steady-state solution, then these can be substituted into the conditions of term dominance for the corresponding Case, thus defining the boundaries for its region of validity (1). The boundary conditions for Case 2, displayed in logarithmic coordinates, are
print 'Boundaries:'
case2.boundaries_log
In general, the boundary conditions are linear in logarithmic space and can be used to formulate linear programming problems. The first phase of a linear programming problem is to determine if the constraints (boundaries) define a feasible region in parameter space (11). This can be tested very efficiently and independent of specific parameter values. If the boundaries define a feasible region, then the case represents a phenotype of the system; otherwise, the case is not a valid phenotype. Case 2 is tested for validity as shown next,
print 'Case is valid? ' + str(case2.is_valid())
In this example we have automatically determined that this case is not a valid phenotype of the system because the boundaries are not simultaneously satisfied anywhere in parameter space. This is easily shown for this simple example by inspecting the boundary conditions, as discussed in Online Methods.
The complete set of potentially valid phenotypes can be analyzed to identify cases that are valid somewhere in parameter space, and the complete set of valid phenotypes defines the phenotypic repertoire of the system (1). Validity of the phenotypes can be determined automatically and independently; thus, they are computed concurrently using parallel computing methods. We calculate the complete phenotypic repertoire of the example system, as shown below.
# We add bounds on the values of the rho parameters to ensure conserved mode of regulation e.g. activation vs repression
phenotypic_repertoire = ds.valid_cases(p_bounds={'rho1':[1, 1e20], 'rho2':[1, 1e20]})
The phenotypic repertoire determines the potential behaviors a system can exhibit. Each phenotype can be analyzed for its steady-state solution, signal amplification factors, parameter sensitivities, local stability, etc. Each phenotype is an S-System that can be analyzed for a wealth of properties, using a variety of methods (12). In Table S1 of the Online Methods we show the complete phenotypic repertoire for this example system, together with the signal amplification factors and local stability of each phenotypes at a representative point.
Here, we reproduce the data in this table automatically by inspecting the DesignSpace and Case objects.
print 'Automatically generated reproduction of Table S1 in the Online Methods.\n'
# We print the table headers.
headers = ['Case #', 'Signature', 'Validity', 'Log gain', 'Stability']
print ' '.join([column.center(10) for column in headers])
# We iterate over all the potentially valid phenotypes
for case_number in range(1, ds.number_of_cases+1):
case = ds(case_number)
if case_number not in phenotypic_repertoire:
# If the case is not part of the phenotypic repertoire we only print '-'
row = [str(case_number), case.signature, '-', '-', '-']
else:
ssys = case.ssystem
# We calculate the signal amplification factor for X2 as a function of X3 (log gain).
amplification_fx = ssys.log_gain('X2', 'X3')
p = case.valid_parameter_set()
# We calculate the number of eigenvalues with positive real parts.
positive_real_eigenvalues = case.positive_roots(p)
if positive_real_eigenvalues == 0:
stability = 'S'
elif positive_real_eigenvalues == 1:
stability = 'U'
# Add these to a list to be printed
row = [str(case_number),
case.signature, '+',
str(amplification_fx).rjust(4),
stability]
print ' '.join([column.center(10) for column in row])
The analysis of the model up until this point has been completely determined by the archiectual components of the model, which includes the connectivity, the rate law equations, the kinetic orders and the non-negativity of \(\log\rho_1\) and \(\log\rho_2\). The analysis we have shown thus far has been enitirely based on queries on the various objects that perform all the necessary calculations internally. This automated approach has given us a wealth of information at a global level that includes the available responses to environmental stimulus and the stability of different phenotypes.
Automatic Identification of Phenotype-Specific Sets of Parameter Values
In this subsection we show the methods for automatically determining sets of parameter values for phenotypes of a system. These parameter sets can be used as test values to get an idea of the phenotypic context, and the values may be refined by bounding them within specific domains of design space. To show this component of our method, we begin by determining a set of parameter values for the phenotype identified as Case 5. We will then explore the context of this phenotype graphically using a 2-D design space plot, centered on the values automatically determined.
First, we acquire the Case object for Case 5 and we display its case signature and its S-System equations.
case5=ds(5)
print 'Case number: ' + str(case5.case_number)
print 'Case signature ' + '[' + case5.signature + ']\n'
print 'Equations: '
case5.equations
Once the Case object has been acquired, we automatically determine a set of parameter values within the region of validity, as shown below.
parameter_values = case5.valid_parameter_set()
for parameter_name in parameter_values:
print parameter_name + ' = ' + str('%.4f' % parameter_values[parameter_name])
Using these parameter values as a reference set, we define a 2-D slice of design space, where two slice parameters are allowed to change freely. In this example, we chose \(X_3\) and \(K_1\) as the slice parameters and for each we define a range for our plots. The range we have selected is 6 decades for each parameter, centered on the automatically determined values.
x_parameter = 'X3'
y_parameter = 'K1'
range_x = [parameter_values['X3']*1e-3, parameter_values['X3']*1e3]
range_y = [parameter_values['K1']*1e-3, parameter_values['K1']*1e3]
We can automatically determine the phenotypes that are valid within some \(n\)-D slice of design space. To achieve this, we create a dictionary of parameter values and ranges that defines the particular slice. For example, we define the following dictionary for a 2-D slice with \(X_3\) and \(K_1\) slice parameters. The non-slice parameters are held fixed at their reference values.
# We create a new dictionary that using the automatically determined values for the parameters.
slice_dict = dict(parameter_values)
# Set the x and y parameters of the dictionary to be ranges specified in the previous block.
slice_dict[x_parameter] = range_x
slice_dict[y_parameter] = range_y
The slice_dict dictionary is passed to the DesignSpace object to determine the valid cases within this slice.
phenotypes_in_slice = ds.valid_cases(p_bounds=slice_dict)
# We display the valid cases as a table with Case # and Signature
print 'Cases in slice:'
print ' Case # \tSignature'
# We iterate over all the phenotypes valid in the slice
for i in phenotypes_in_slice:
case = ds(i)
print ' '+str(case.case_number)+' \t'+case.signature
This can be shown graphically by plotting the valid regions of each case by using the plotting utilities included in the dspace.plotutils module. We create a new figure using matplotlib.pyplot.figure and create an axes object that serves as the drawing canvas. These object are passed as arguments to the draw_2D_slice method of the DesignSpace object, together with the necessary information to define the 2-D slice,
# These colors are used in Figure S2 of Online Methods, and are added here for consistency
color_dict = {'1': (1.0, 0.0, 0.16, 1.0),
'5': (1.0, 0.43137254901960792, 0.0, 1.0),
'15': (1., 1., 0., 1.),
'16': (0.0, 1.0, 0.84334809192494242, 1.0)
}
fig = pyplot.figure()
ax = pyplot.gca()
cdict=ds.draw_2D_slice(ax, parameter_values,
x_parameter, y_parameter,
range_x, range_y,
color_dict=color_dict,
colorbar='single')
Each region corresponds to a qualitatively distinct phenotype of the system. These regions are obtained through vertex enumeration by virtue of the linear boundaries in logarithmic space.
For each of the cases in this slice we can acquire a set of parameter values. To enforce that the automatically determined values for the parameters be within this slice, we impose bounds on the values the parameters can assume. Furthermore, we can identify parameter sets that are within the interior of the feasible region.
Sets of parameter values in the interior of the feasible region can be automatically found by specifying a desired fold-difference from the nearest boundaries in each orthogonal direction. As this fold-difference increases, the point we obtain tends towards the center of the region. In the next series of steps we generate an array of values indicating desired fold-differences and automatically determine an array of parameter values with increasing fold-differences from the boundaries.
# We specify an array of values that fold-differences from the boundaries.
distance=np.logspace(np.log10(2), 3, 10)
print distance
For each valid case in this slice, we determine multiple parameter sets that converge at the centroid of the feasible region. Each set of parameter values is determined by requiring a defined fold-difference from the nearest boundaries (using the values shown above). Note, if the nearest boundaries are less than the required fold-difference from the starting point in both directions, then the geometric mean of the distances between the boundaries is chosen.
# Define a dictionary that will hold the set of x and y values
# for each of the phenotypes in the slice.
point_sets = dict()
# We iterate over the set of phenotypes in the slice
for case_number in phenotypes_in_slice:
# Acquire the case objects
case=ds(case_number)
# Acquire a set of parameter values in this slice (the starting set)
case_parameters=case.valid_parameter_set(p_bounds=slice_dict)
# Store the value of the x and y parameters in a list
points=[np.log10([case_parameters[x_parameter],
case_parameters[y_parameter]])]
for k in distance:
# For each value in the distance array we calculate a new set of parameter values
case_parameters2=case.valid_interior_parameter_set(p_bounds=slice_dict,
distance=k)
# We store the values for the x and y parameters
points.append(np.log10([case_parameters2[x_parameter],
case_parameters2[y_parameter]]))
point_sets[case_number] = points
Each parameter set that is acquired using this strategy corresponds to a single point in this 2-D space. These points can be overlayed on the 2-D design space plot shown previously. The sequence of points generated as the fold-difference from the boundaries increases creates a trajectory from the edge of the feasible region towards the center. The final parameter set, with the greatest defined distance from the boundaries, is represented by a square in the Figure below.
fig = pyplot.figure()
ax = pyplot.gca()
cdict=ds.draw_2D_slice(ax, parameter_values,
x_parameter, y_parameter,
range_x, range_y,
color_dict=color_dict,
colorbar='single')
for case_number in phenotypes_in_slice:
points = point_sets[case_number]
ax.plot([i[0] for i in points], [i[1] for i in points], 'ko-',
mfc='w', mec='k', ms=2.5)
ax.plot(points[-1][0], points[-1][1], 's',
mfc='w', mec='k', ms=4)
This strategy for identifying parameter values is a powerful tool to explore the phenotypic landscape of biochemical models. Parameter sets can be efficiently acquired for phenotypes with very small feasible regions in high dimensional spaces that may easily be overlooked if we rely on sampling of parameter space. Through this method we can find parameter values for all the phenotypes and these parameter sets can serve to characterize the system's behavior using conventional methods.
Automatic Identification of Parameter Sets for Co-Localized Phenotypes
In the previous subsection we focused on identifying parameter values for single phenotypes. However, using the design space approach we can find parameter values such that multiple phenotypes are found together in some slice of design space. In this section we show the methods for identifying slices where phenotypes of the system are co-localized, meaning that they are simultaneously realized within some slice of design space.
First, we focus on the strategy for identifying the greatest number of phenotypes that can be co-localized within a slice of design space. In the second subsection, we focus on identifying slices of design space with specific ensembles of phenotypes.
Maximizing the Number of Phenotypes in a Slice of Design Space
The first problem we show is identifying the maximum number of phenotypes for a set of slice variables. To achieve this, we calculate the phenotypic repertoire of the system and then determine which combinations of these yields the largest number of co-localized phenotypes. The strategy can be optimized to reduce the number of combinations tested, although the problem remains combinatorial and the number of combinations to check can be astronomical for larger systems. Nonetheless, it has useful applications for smaller systems, as shown in the following example.
We select \(K_1\) and \(K_2\) as the slice parameters and check for the maximum number of co-localized phenotypes given the entire phenotypic repertoire.
co_localize_phenotypes = ds.maximum_co_localized_cases(['K1', 'K2'],
phenotypic_repertoire)
print 'Maximum number of co-localized phenotypes: ' + str(len(co_localize_phenotypes[0]))
If we compare this result with those in Table S3, we see that there is a slice such that the entire phenotypic repertoire is simultaneously realized (co-localized), and a set of parameters can be determined automatically.
sets_of_parameters = ds.co_localize_cases(co_localize_phenotypes[0], ['K2', 'K1'])
This method returns a dictionary that contains pairs of case numbers and a set of parameter values that places it in the slice for co-localization. The sets of parmeter values for this example are shown in the table below. Note that the only values that are different between the cases are those of the slice parameters, \(K_1\) and \(K_2\).
paramaters_case13 = sets_of_parameters['13']
print '\nTable of values for K1 and K2'
headers = ['Parameter'] + ['Case '+str(i) for i in co_localize_phenotypes[0]]
print ' '+' '.join([i.center(10) for i in headers])
for name in ['K1', 'K2']+[i for i in paramaters_case13.keys() if i not in ['K1', 'K2']]:
row = [name]
for case_number in co_localize_phenotypes[0]:
# Format to include 2 significant digits
row += ['%.5f'% sets_of_parameters[str(case_number)][name]]
print ' '+' '.join([i.center(10) for i in row])
We represent this slice graphically by plotting the design space using these parameters as our reference. We select a case, here we use Case 13, and we select bounds on the slice parameters and apply the plotting methods.
x_parameter = 'K1'
y_parameter = 'K2'
range_x = [paramaters_case13['K1']*1e-2, paramaters_case13['K1']*1e2]
range_y = [paramaters_case13['K2']*1e-2, paramaters_case13['K2']*1e2]
fig = pyplot.figure()
ax = pyplot.gca()
cdict=ds.draw_2D_slice(ax, paramaters_case13,
x_parameter, y_parameter,
range_x, range_y,
intersections=[1, 3],
color_dict=color_dict,
colorbar='single')
xticks=ax.set_xticks([-2, -1., 0, 1, 2])
yticks=ax.set_yticks([-2, -1., 0, 1, 2])
Inspection of this figure shows that there are total of 9 regions, 6 of which correspond to regions with a unique Case, and 3 that correspond to regions of overlap of multiple Cases. The regions of overlap are the signature for regions of multistability. Note that 3 phenotypes do not appear by themselves and appear strictly within these regions of overlap and correspond to the Cases that are the dynamically unstable. The ability to identify the slice with the most phenotypes is often useful to understand the phenotypic potential of a system, and to identify parameter values that result in interesting phenotypic landscapes.
Ensembles of Desireable Phenotypes
Case co-localization can also be used to identify slices that exhibit particular ensembles of phenotypes. In the following blocks we identify a slice with a specific progression of behaviors as a slice parameter, in this case \(X_3\), increases.
For example, assume that we are interested in finding parameter values that yield the following sequence of behaviors:
At low levels of \(X_3\) there is a basal level of \(X_2\); as \(X_3\) increases, the system enters a regulatable regime where changes in \(X_2\) are directly coupled to changes in \(X_3\) until it enters a bistable regime that creates a hysteretic loop; when \(X_3\) reaches a threshold, it switches to a higher state that remains regulatable (\(X_2\) directly coupled to \(X_3\)); finally, the effect of \(X_3\) saturates at a maximal level of \(X_2\) (see Online Methods for detailed description).
To achive this sequence of phenotypes, we (a) select the slice parameter, (b) select the cases with the appropriate behaviors from the phenotypic repertoire to define the ensemble to be co-localized, (c) define the constraints that specifies a particular sequence of behaviors and (d) automatically determine the parameter values that give this progression of phenotypes.
In this example, we choose the phenotypes that have the correct signal amplification factors and stability of the steady states: Cases 1, 5, 13, 15, and 16.
slice_parameter = 'X3'
ensemble_cases = [1, 5, 13, 15, 16]
Next, we define constraints between the slice variables for each of the Cases in the ensemble. Note: The notation of these constraints differs from the notation in the Online Methods. The notation here is functional, whereas the notation in the article is mathematical. Here, we define special auxiliary variables that relate to the slice parameters for each of the Case objct using the following notation,
'$<slice parameter name>_<case index in list of ensemble cases>'.
Therefore, \(s_1\) in Eq S31 in the Main Text is input here as $X3_0 in arithmetic space. The constraints must be positive powerlaws on both sides of the inequality, defined again using a string representation.
# We define the sequence of phenotypes corresponding to Eqs. S31-S34 in Online Methods
sequence_constraints = ['$X3_0 < $X3_1',
'$X3_1 < $X3_2',
'$X3_2 < $X3_3',
'$X3_3 < $X3_4']
# We get a set of parameter values for each case using these constraints
sets_of_parameters = ds.co_localize_cases(ensemble_cases,
slice_parameter,
constraints=sequence_constraints)
The slice is represented graphically by plotting the logarithm of the concentration of \(X_3\) on the \(x\)-axis and the logarithm of the steady-state concentration of \(X_2\) on the \(y\)-axis. The stability of the steady states is shown by the line style, where a solid line represents a stable steady state and a dash-dot line represents an unstable steady states. Colored dots are added to the plot showing the automatically determined values for each of the phenotypes.
paramaters_case13 = sets_of_parameters['13']
range_x = [paramaters_case13['X3']*10**-2.1, paramaters_case13['X3']*10**2.1]
fig = pyplot.figure()
ax1 = fig.add_axes([0.1, 0.1, 0.7, 0.8])
ds.draw_1D_positive_roots(ax1, 'log(X2)', paramaters_case13,
slice_parameter, range_x,
line_dict={1:{'ls':'-.', 'c':'k'},
0:{'ls':'-', 'c':'k'}})
ax1.set_ylim([-6, 6])
ax1.set_xlim(np.log10(range_x))
ax1.set_ylabel('$\log_{10} X_2$')
ax1.set_xlabel('$\log_{10} X_3$')
xticks=ax1.set_xticks([np.log10(paramaters_case13['X3'])+i for i in [-2, -1, 0, 1, 2]])
for i in sets_of_parameters:
pvals = sets_of_parameters[i]
x = sets_of_parameters[i][slice_parameter]
case = ds(i)
y = case.ssystem.steady_state_function('log(X2)', sets_of_parameters[i])
ax1.plot(np.log10(x), y, '.', markerfacecolor=cdict[i], markeredgecolor='k')
ax2 = fig.add_axes([0.85, 0.1, 0.1, 0.8])
ds.draw_region_colorbar(ax2, {i:cdict[i] for i in sets_of_parameters})
Another example of case co-localization is shown in the next few block, where we identify a 1-D slice such that all the unstable phenotypes are co-localized. Again, we select \(X_3\) as the slice parameter.
To achive this ensemble of phenotypes, we (a) select the slice parameter, (b) select the cases in the ensemble to be co-localized and (c) automatically determine the parameter values that give this ensemble. Note, here we do not require a specific sequence of phenotypes and hence do not specify additional constraints.
slice_parameter = 'X3'
ensemble_cases = [9, 13, 14]
sets_of_parameters = ds.co_localize_cases(ensemble_cases, slice_parameter)
The slice is represented graphically as in the previous plot. Here, the plot shows all unstable phenotypes in a 1-D slice, and each unstable phenotype is sandwiched between two stable phenotypes, which is indicative of bi-stability throughout the entire slice.
paramaters_case13 = sets_of_parameters['13']
range_x = [paramaters_case13['X3']*10**-2.1, paramaters_case13['X3']*10**2.1]
fig = pyplot.figure()
ax1 = fig.add_axes([0.1, 0.1, 0.7, 0.8])
ds.draw_1D_positive_roots(ax1, 'log(X2)', paramaters_case13,
slice_parameter, range_x,
line_dict={1:{'ls':'-.', 'c':'k'},
0:{'ls':'-', 'c':'k'}}
)
ax1.set_ylim([-6, 3])
ax1.set_ylabel('$\log_{10} X_2$')
ax1.set_xlabel('$\log_{10} X_3$')
yticks=ax1.set_yticks([-6, -3, 0, 3])
xticks=ax1.set_xticks([-2, -1, 0, 1, 2])
for i in sets_of_parameters:
pvals = sets_of_parameters[i]
case = ds(i)
x = sets_of_parameters[i][slice_parameter]
y = case.ssystem.steady_state_function('log(X2)', sets_of_parameters[i])
ax1.plot(np.log10(x), y, '.', markerfacecolor=cdict[i], markeredgecolor='k')
ax2 = fig.add_axes([0.85, 0.1, 0.1, 0.8])
ds.draw_region_colorbar(ax2, {i:cdict[i] for i in sets_of_parameters})
Together, these examples show the ability of the design space framework to quickly and efficiently find specific phenotypic landscapes in high dimensional spaces. This method has potential applications in the analysis of both natural and synthetic circuits alike in which specific responses to environmental stimuli are observed or desired.
Application to a General Class of Oscillator Designs
In the previous section we emphasized the methods using a simpler example. Here, we focus on systems that have more variables and more complex regulatory mechanisms. We begin with a two-gene circuit that encodes a relaxation oscillator that has been shown to display oscillations (13, 3). Then, we extend the analysis to a general class of two-gene circuits with several designs that we have previously analyzed and compared with regard to the global robustness of their oscillatory phenotypes (10). The analysis we show here does not assume or specify values for any of the parameters.
Application to a Relaxation Oscillator
We begin with the analysis of the relaxation oscillator. The analysis of a new system requires that we construct a new instance of the DesignSpace object. The input to the dspace package must be in GMA form, thus we use Eqs. S39-S44 in Online Methods. The recast equations are input into the dspace package as strings, which we use to instantiate the Equations and consequently the DesignSpace objects. In this example, we explicitly add the constraints to the design space object where \(\rho_1\) and \(\rho_3\) must be greater than one, to ensure the qualitative nature of transcriptional control is not altered.
equation_strings = ['X1. = a1*rho1^-1*X5^-1 \
+ a1*X2^2*K2A^-2*X5^-1 \
+ a1*rho1^-1*X4^2*K4^-2*X5^-1 - b1*X1',
'X2. = a2*X1 - b2*X2',
'X3. = a3*rho3^-1*X6^-1 + a3*X2^2*K2R^-2*X6^-1 - b3*X3',
'X4. = a4*X3 - b4*X4',
'X5 = 1 + X2^2*K2A^-2 + X4^2*K4^-2',
'X6 = 1 + X2^2*K2R^-2'
]
# Again we define the simplified to latex dictionary for the model
symbols = {'X1':'X_1',
'X2':'X_2',
'X3':'X_3',
'X4':'X_4',
'X5':'X_5',
'X6':'X_6',
'a1':r'\alpha_1',
'a2':r'\alpha_2',
'a3':r'\alpha_3',
'a4':r'\alpha_4',
'b1':r'\beta_1',
'b2':r'\beta_2',
'b3':r'\beta_3',
'b4':r'\beta_4',
'rho1':r'\rho_1',
'rho3':r'\rho_3',
'K2R':'K_{2R}',
'K2A':'K_{2A}',
'K4':'K_4'}
equations_object = dspace.Equations(equation_strings, latex_symbols=symbols)
ds = dspace.DesignSpace(equations_object, constraints=['rho1 > 1', 'rho3 > 1'])
The recast equations for this system are shown below,
ds.equations
we automatically enumerate the model phenotypes and identify the phenotypic repertoire of the system,
phenotypic_repertoire = ds.valid_cases()
with the phenotypic repertoire identified, we automatically determine a set of parameter values and determine if the phenotype is stable, unstable, or has the potential to oscillate. We achieved this using a simple strategy: we (1) automatically determine a set of parameter values within the interior of this phenotype and (2) calculate the number of eigenvalues with a positive real part using the Routh criteria (14). Using these two steps, we reproduce the data shown in Table 1 of the Main Text,
print 'Automatically generated reproduction of Table S1.\n'
# We define the table column headers in a list.
headers = ['Case #', 'Signature', '# of pos. real eigen.']
widths = [6, 12, 23]
print ' '+' '.join([headers[i].center(widths[i]) for i in range(len(headers))])
for case_number in phenotypic_repertoire:
case = ds(case_number)
p = case.valid_parameter_set()
# We calculate the number of eigenvalues with positive real parts.
positive_real_eigenvalues = case.positive_roots(p)
row = [str(case_number),
case.signature,
str(positive_real_eigenvalues)]
print ' '+' '.join([row[i].center(widths[i]) for i in range(len(row))])
From the table above we have found one phenotype that has the potential to oscillate (two eigenvalues with positive real part).
In the following block we explore the phenotype with oscillatory potential, Case 23, and its context by looking at the neighboring phenotypes. We create the Case object corresponding to Case 23 and automatically determine a set of parameter values that is 3 decades away from the nearest boundaries.
case23 = ds(23)
pvals = case23.valid_interior_parameter_set(distance=1e3)
for parameter_name in pvals:
print parameter_name + ' = ' + str(pvals[parameter_name])
We visualize this phenotype and its context with a 2-D design space plot. We specify the slice parameters, \(\beta_2\) and \(\beta_4\), and a range for each to plot, in this example 2 decades centered on the automatcally determined parameter values for Case 23 (●).
x_parameter = 'b2'
y_parameter = 'b4'
range_x = [pvals[x_parameter]*1e-2, pvals[x_parameter]*1e2]
range_y = [pvals[y_parameter]*1e-2, pvals[y_parameter]*1e2]
fig = pyplot.figure()
fig.set_size_inches(6, 3)
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
cdict=ds.draw_2D_slice(ax, pvals,
x_parameter, y_parameter,
range_x, range_y, intersections=[1, 3, 5])
axticks=ax.set_xticks([-2, -1, 0, 1, 2])
axticks=ax.set_yticks([-2, -1, 0, 1, 2])
_=ax.plot(np.log10(pvals[x_parameter]), np.log10(pvals[y_parameter]), 'k.')
We can determine the stability of the fixed points within this slice, as determined by the number of eigen-values with a positive-real part, and represent this visually in the following plot,
# We define a color dictionary for the stability plots.
stability_colors = {'0':(0.0, 0.8471, 0.9843, 1.0),
'2':(1., 1., 0., 1.),
'1':(1., 0., 0., 1.),
'0,1':(1., 0., 0., 1.),
'0,1,2':(1., 0.5, 0., 1.)}
fig = pyplot.figure()
fig.set_size_inches(3, 1.8)
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
cdict=ds.draw_2D_positive_roots(ax, pvals,
x_parameter, y_parameter,
range_x, range_y,
color_dict=stability_colors,
resolution=300)
_=ax.plot(np.log10(pvals[x_parameter]), np.log10(pvals[y_parameter]), 'k.')
_=ax.set_xticks([-2, -1, 0, 1, 2])
_=ax.set_yticks([-2, -1, 0, 1, 2])
This plot shows regions of monostability (blue), oscillatory potential (yellow), multistability (red) and overlaps between stable, exponentially unstable and oscillatory unstable phenotypes (orange). The orange regions have the potential to generate global bifurcations known as Saddle Node into Limit Cycle bifurcations (SNIC) that result in a transition from a stable steady state to full-blown oscillations (10).
The potential for oscillation can be tested by simulating the original equations. Here, we simulate the dynamics of the full system using a thin python interface for the Sundials CVODE solver (15). We import the python interface to the CVODE solver into the python environment. This interface only accepts systems of differential equations; in order to simulate the system we define the original equations using a string representation.
from integrate import integrate
original_system =['X1. = a1*(rho1^-1 + X2^2*K2A^-2 + rho1^-1*X4^2*K4^-2) \
/(1 + X2^2*K2A^-2 + X4^2*K4^-2) - b1*X1',
'X2. = a2*X1 - b2*X2',
'X3. = a3*(rho3^-1 + X2^2*K2R^-2)/(1 + X2^2*K2R^-2) - b3*X3',
'X4. = a4*X3 - b4*X4'
]
Once the original system is defined in this way, we process it for input into the integrate module. This module accepts a dictionary, where the keys are the names of the time-differentiated variables and the values are the string representation of the right-hand side of the differential equation. Thus, we split the string representation of the original system and retain the right-hand side of the equations.
equation_rhs = [equation.split(' = ')[1] for equation in original_system]
The integrate module uses python to evaluate the right hand side of the equations. Python uses '**' as the power operator, which must replace all occurences of '^' in the original equations.
equation_rhs = [equation.replace('^','**') for equation in equation_rhs]
Using these equations we define the dictionary of variable-equation pairs,
equation_dictionary = {ds.dependent_variables[i]:equation_rhs[i] for i in range(len(equation_rhs))}
We choose the steady-state values at the automatically determined parameter values obtained from the design space approximation as the initial conditions for the simulation.
initial_conditions = case23.steady_state(pvals)
# We remove the values for the auxiliary variables (not defined for the original systems).
_=initial_conditions.pop('X5')
_=initial_conditions.pop('X6')
These are approximate values but they should be near the actual steady-state solution and the instability will repel the system from the actual steady-state solution towards a stable limit cycle, as shown below,
t, x = integrate(equation_dictionary, initial_conditions, pvals, 0, 200, 0)
fig = pyplot.figure()
fig.set_size_inches(3., 1)
ax = fig.add_axes([0.1, 0.1, 0.7, 0.7])
ax.plot(t, np.log10(x['X4']), 'k')
ax.set_ylim(-2, 0)
_=ax.set_yticks([-2, -1, 0])
_=ax.set_xlabel('$t$')
_=ax.set_ylabel('$\log_{10}[X_4]$')
The simulation shows that the potential for oscillation results in true oscillations. The oscillations are centered around the value of the initial conditions that should be near the steady-state solution of the system. If we keep the \(\beta_4\) parameter fixed at \(\beta_4 = 1\), and allow \(\beta_2\) to vary, we see a series of dynamic behaviors. According to the design space approximation, as we increase \(\beta_2\) from \(\beta_2 = 1\times10^{-2}\), the system goes from being stable, then bistable (hysteretic), then oscillatory, then overlap between unstable, stable and oscillatory, followed by stable. The induction curve for the original system can be shown using a simulations using a strategy involving numerical continuation.
We simulate the system at \(\beta_2 = \beta_{2,0}\) until it reaches an attractor (steady state or limit cycle) and record the min and max values. The final values of the simulation are stored as the initial conditions for the next simulation with a small increase in the value of the \(\beta_2\) parameter, given by \(\Delta\beta_2\), and again we record the min and max values of the attractor. This is repeated iteratively such that the value for \(\beta_2\) at the \(i\)-th iteration is \(\beta_{2,i} = \beta_{2,i-1}+\Delta\beta_2\), where \(i = {1, 2, 3, ..., n}\). We repeat this process starting with \(\beta_2 = \beta_{2,n}\) and reverse the direction so \(\beta_{2,i} = \beta_{2,i+1} - \Delta\beta_2\) until \(\beta_2 = \beta_{2,0}\).
pvals['b2'] = 1e-2
pvals['b4'] = 1
case = ds(ds.valid_cases(p_bounds=pvals)[0])
initial_conditions = case.steady_state(pvals)
# We remove the values for the auxiliary variables (not defined for the original systems).
_=initial_conditions.pop('X5')
_=initial_conditions.pop('X6')
x_vals1 = np.logspace(-2, 1, 200)
y_vals1 = dict(min_val = [], max_val = [])
x_vals2 = np.logspace(1, -2, 200)
y_vals2 = dict(min_val = [], max_val = [])
for x in x_vals1:
pvals['b2'] = x
t, y = integrate(equation_dictionary, initial_conditions, pvals, 0, 2000, 0)
initial_conditions = {i:y[i][-1] for i in y}
t, y = integrate(equation_dictionary, initial_conditions, pvals, 0, 1000, 0)
y = y['X2']
y_vals1['min_val'].append(np.min(y))
y_vals1['max_val'].append(np.max(y))
for x in x_vals2:
pvals['b2'] = x
t, y = integrate(equation_dictionary, initial_conditions, pvals, 0, 2000, 0)
initial_conditions = {i:y[i][-1] for i in y}
t, y = integrate(equation_dictionary, initial_conditions, pvals, 0, 1000, 0)
y = y['X2']
y_vals2['min_val'].append(np.min(y))
y_vals2['max_val'].append(np.max(y))
The results of the forward and reverse continuations are shown below. The forward continuation is shown in red, the backward continuation is shown in blue. An overlap of the forward and reverse continuations is shown in purple. The approximated bifurcation using the design space method is represented by the gray line: a solid line represents a stable steady state, dashed line represents an oscillatory steady state and the dotted line represents an unstable steady state.
fig = pyplot.figure()
fig.set_size_inches(3., 1.3)
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
ds.draw_1D_positive_roots(ax, 'log(X2)', pvals,
'b2',
[1e-2, 1e1],
resolution=200,
line_dict={2:{'ls':'--', 'c':'gray', 'lw':0.7},
1:{'ls':':', 'c':'gray', 'lw':0.7},
0:{'ls':'-', 'c':'gray', 'lw':0.7}}
)
ax.fill_between(np.log10(x_vals1), np.log10(y_vals1['min_val']), y2=np.log10(y_vals1['max_val']), color='r', alpha=0.5)
ax.fill_between(np.log10(x_vals2), np.log10(y_vals2['min_val']), y2=np.log10(y_vals2['max_val']), color='b', alpha=0.5)
ax.set_xlabel(r'$\log_{10}(\beta_2)$')
ax.set_ylabel(r'$\log_{10}[X_2]$')
ax.set_ylim([-4, 4])
ax.set_yticks([-4, -2, 0, 2, 4])
_=ax.set_xticks([-2, -1, 0, 1])
The numerical continuation shows regions of bistability (forward and reverse continuations do not overlap) that result in a hysteretic loop. In addition, we see a region of oscillation given by the envelope of the oscillation (shown by the purple area where the min and max do not overlap) that begins with small amplitudes (super-critical hopf bifurcation) and abruptly collapses to a stable node (SNIC bifurcation). In addition, the design space approximation correctly predicts the existance of stable, unstable and oscillatory phenotypes but over-estimates the regions of hysteresis and oscillation (3).
Maximizing Number of Phenotypes in a 2-D Slice
In the final example for this relaxation oscillator, we identify the maximum number of phenotypes that can be co-localized on a 2-D slice of design space, where the slice parameters are \(\beta_2\) and \(\beta_4\).
co_localize_phenotypes = ds.maximum_co_localized_cases(['b2', 'b4'],
ds.valid_cases())
for i in co_localize_phenotypes:
print 'Slice with ' + str(len(i)) + ' cases. Cases in slice: ' + str(i)
Our results indicate that there are three distinct sets of 11 phenotypes that can be co-localized in a 2-D slice, as shown above. The example design space plot for this system, shown in Figure 2 of the Main Text, also consists of 11 phenotypes and involves those in the third set. Thus, we have already shown an example of a slice that maximizes the number of phenotypes. Although this is just one of the three possible sets, examples for the other two sets can be readily determined.
General Class of Two-Gene Oscillator Designs
Here, we extend the analysis shown in the previous subsection to a general class of 16 unique two-gene circuits. As in the previous example, we focus on phenotypes that have the potential to exhibit oscillatory behaviors. This class of circuits can be represented with a single mathematical model that we have previously formulated in an analysis of 7 of these 16 designs (10). The recast equations shown in Eqs. S51-S58 in Online Methods are used to construct the Equations object.
equations_general = ['x1. = B10*(a1/a10)*rho10^0.5*rho1^-pi1*X5^-1\
+ B10*(a1/a10)*rho10^0.5*(K20^g12)*delta1*x2^g12*K2^-g12*X5^-1 \
+ B10*(a1/a10)*rho10^0.5*(K40^g14)*rho1^-1*x4^g14*K4A^-g14*X5^-1\
- b1*x1',
'xA. = BA0*(aA/aA0)*x1 - bA*xA',
'x2. = B20*(bA/bA0)*(n20*a10)*(K20^-1*B10^-1*rho10^-0.5)*xA - b2*x2',
'x3. = B30*(a3/a30)*rho30^0.5*rho3^-pi3*X6^-1\
+ B30*(a3/a30)*rho30^0.5*(K20^g32)*x2^g32*K2R^-g32*X6^-1\
+ B30*(a3/a30)*rho30^0.5*(K40^g34)*delta3*rho3^-1*x4^g34*K4^-g34*X6^-1\
- b3*x3',
'xR. = BR0*(aR/aR0)*x3 - bR*xR',
'x4. = B40*(bR/bR0)*(n40*a30)*(K40^-1*B30^-1*rho30^-0.5)*xR - b4*x4',
'X5 = 1 + (K20^g12)*delta1*x2^g12*K2^-g12 + (K40^g14)*x4^g14*K4A^-g14',
'X6 = 1 + (K20^g32)*x2^g32*K2R^-g32 + (K40^g34)*delta3*x4^g34*K4^-g34']
Again, we define the dictionary relating the simplified notation with the latex notation for each of the parameters and variables in the model.
symbols = {'x1':'x_1', 'xA':'x_A', 'x2':'x_2',
'x3':'x_3', 'xR':'x_R', 'x4':'x_4',
'X5':'X_5', 'X6':'X_6', 'a1':r'\alpha_1',
'a2':r'\alpha_2', 'aA':r'\alpha_A', 'a3':r'\alpha_3',
'aR':r'\alpha_R', 'a4':r'\alpha_4', 'a10':r'(\alpha_1^0)',
'aA0':r'(\alpha_A^0)', 'a20':r'(\alpha_2^0)', 'a30':r'(\alpha_3^0)',
'aR0':r'(\alpha_R^0)', 'a40':r'(\alpha_4^0)', 'B1':r'B_1',
'BA':r'B_A', 'B2':r'B_2', 'B3':r'B_3', 'BR':r'B_R',
'B4':r'B_4', 'bA':r'\beta_A', 'bR':r'\beta_R', 'bA0':r'(B_A^0)',
'bR0':r'(B_R^0)', 'B10':r'(B_1^0)', 'BA0':r'(B_A^0)', 'B20':r'(B_2^0)',
'b1':r'\beta_1', 'b2':r'\beta_2', 'b3': r'\beta_3', 'b4':r'\beta_4',
'B30':r'(B_3^0)', 'BR0':r'(B_R^0)', 'B40':r'(B_4^0)', 'rho1':r'\rho_1',
'rho3':r'\rho_3', 'rho10':r'(\rho_1^0)', 'rho30':r'(\rho_3^0)', 'K2':'K_2',
'K4':'K_4', 'K20':'(K_2^0)', 'K40':'(K_4^0)', 'n20':'n_2^0', 'n40':'n_4^0',
'K2R':'K_{2R}', 'K4A':'K_{4A}', 'delta1':'\delta_1', 'delta3':'\delta_3',
'pi1':'\pi_1', 'pi3':'\pi_3', 'g12' : 'g_{12}','g14' : 'g_{14}',
'g32' : 'g_{32}', 'g34' : 'g_{34}'
}
equations_object = dspace.Equations(equations_general, latex_symbols=symbols)
For convenience, the recast equations for the general class of two-gene circuits are displayed below,
equations_object
Each design is obtained from these equations by choosing a unique combination of values for the \(\delta_1\), \(\pi_1\), \(\delta_3\) and \(\pi_3\) binary variables. A unique combination can be stored in a dictionary that is associated to a specific system, as shown by the following nested dictionaries,
substitutions = {'S.1': {'delta1': '0', 'pi3': '0', 'delta3': '0', 'pi1': '0'},
'S.2': {'delta1': '1', 'pi3': '0', 'delta3': '0', 'pi1': '0'},
'S.3': {'delta1': '0', 'pi3': '0', 'delta3': '1', 'pi1': '0'},
'S.4': {'delta1': '1', 'pi3': '0', 'delta3': '1', 'pi1': '0'},
'S.5': {'delta1': '0', 'pi3': '0', 'delta3': '0', 'pi1': '1'},
'S.6': {'delta1': '1', 'pi3': '0', 'delta3': '0', 'pi1': '1'},
'S.7': {'delta1': '0', 'pi3': '0', 'delta3': '1', 'pi1': '1'},
'S.8': {'delta1': '1', 'pi3': '0', 'delta3': '1', 'pi1': '1'},
'S.9': {'delta1': '0', 'pi3': '1', 'delta3': '0', 'pi1': '0'},
'S.10':{'delta1': '1', 'pi3': '1', 'delta3': '0', 'pi1': '0'},
'S.11':{'delta1': '0', 'pi3': '1', 'delta3': '1', 'pi1': '0'},
'S.12':{'delta1': '1', 'pi3': '1', 'delta3': '1', 'pi1': '0'},
'S.13':{'delta1': '0', 'pi3': '1', 'delta3': '0', 'pi1': '1'},
'S.14':{'delta1': '1', 'pi3': '1', 'delta3': '0', 'pi1': '1'},
'S.15':{'delta1': '0', 'pi3': '1', 'delta3': '1', 'pi1': '1'},
'S.16':{'delta1': '1', 'pi3': '1', 'delta3': '1', 'pi1': '1'}
}
Using these dictionaries to substitute the values for the binary variables, we create the DesignSpace objects for each of the design in this general class. The DesignSpace objects are stored in a dictionary and are accesible by their system key (e.g 'S.1').
Note: Constraints are imposed on the \(\beta\) and \(B\) parameter by virtue of their definition in the model. The \(\beta\) parameters must be greater than the growth rate, \(\mu = 0.6931\), and the parameter \(B\) must be greater than \(\beta + \mu\); therefore, \(\beta > 0.6931\) and \(B > 2\times0.6931\).
systems = dict()
for i in range(16):
# Construct the keys: S.1, S.2, S.3, ..., S.16
key = 'S.'+str(i+1)
# Get the substitution dictionary for each system
subs = dict(substitutions[key])
# Set the value of the exponents equal to 2
subs.update(g12=2, g14=2, g32=2, g34=2)
# Construct the design space objects
ds = dspace.DesignSpace(equations_object,
parameter_dict=subs,
name=key,
constraints=['b1 > 0.6931471805599453',
'b2 > 0.6931471805599453',
'bA > 0.6931471805599453',
'b3 > 0.6931471805599453',
'bR > 0.6931471805599453',
'b4 > 0.6931471805599453',
'B10 > 2*0.6931471805599453',
'B20 > 2*0.6931471805599453',
'BA0 > 2*0.6931471805599453',
'B30 > 2*0.6931471805599453',
'BR0 > 2*0.6931471805599453',
'B40 > 2*0.6931471805599453',
'rho1 > 1',
'rho3 > 1'],
latex_symbols=symbols)
# Append them to the list
systems[key] = ds
In the next blocks we generate the data to reproduce Table 2 in the Main Text.
To reproduce this table we (a) calculate the phenotypic fraction (the number of valid phenotypes over potentially valid phenotypes) and (b) identify the phenotypes that have the potential for oscillation. A phenotype that has oscillatory potential may not have the potential to oscillate in its entire phenotypic region. This makes it a challenge to find the oscillatory phenotypes. Here, we perform a simple sampling that is focused on parameter value calculated automatically within the interior of the phenotypic regions. First, we automatically determining a parameter set in the interior of the valid region. Second, we specify a 4-D slice with \(\beta_2\), \(\beta_4\), \(\beta_A\) and \(\beta_R\) as the slice parameters, extending 4 decades above and below the automatically determined values. Third, we define a uniform grid of \(3\times3\times3\times3 = 81\) points. Fourth, we test if the point is within the valid region. Fifth, we determine if the phenotype has two eigenvalues with positive real part and if so, we add the phenotype to the set of oscillating phenotypes.
The slice parameters were chosen because these contribute to the effective delay within the circuit, which aside from cooperativity (an architectural feature of the model) can influence if a system oscillates.
slice_parameters = ['b2', 'b4', 'bA', 'bR']
slice_range = [1e-4, 1e4]
number_samples = 3
number_oscillating = {}
# For each system
for system_index in range(len(systems)):
key = 'S.'+str(system_index+1)
ds = systems[key]
oscillating_cases = set()
number_oscillating[ds.name] = oscillating_cases
# For each phenotype in the repertoire:
for case_number in ds.valid_cases():
# First step
case = ds(case_number)
# Second step
pvals = case.valid_parameter_set()
centered = [[pvals[i]*slice_range[0], pvals[i]*slice_range[1]] for i in slice_parameters]
# Third step
X = [np.linspace(log10(i[0]), log10(i[1]), number_samples) for i in centered]
for x1 in X[0]:
pvals[slice_parameters[0]] = 10**x1
for x2 in X[1]:
pvals[slice_parameters[1]] = 10**x2
for x3 in X[2]:
pvals[slice_parameters[2]] = 10**x3
for x4 in X[3]:
pvals[slice_parameters[3]] = 10**x4
if case.is_valid(p_bounds=pvals):
continue
roots = case.positive_roots(pvals)
if roots == 2:
oscillating_cases.add(case_number)
break
if case_number in oscillating_cases:
break
The results are shown in the following reproduction of Table 2 of the Main Text,
print 'Automatically generated reproduction of Table S2.\n'
# We define the table column headers in a list.
headers = ['System', ' Mode of Control ', 'Phen. Fraction', '# Osc. Phenotypes']
# widths = [6, 15, 14, 17]
print ' '+' '.join([headers[i] for i in range(len(headers))])
for i in range(16):
ds = systems['S.'+str(i+1)]
subst = substitutions[ds.name]
row = [ds.name.center(len(headers[0]))]
row += [('π1=' + str(subst['pi1'])\
+ ', δ1=' + str(subst['delta1'])\
+ ', π3=' + str(subst['pi3'])\
+ ', δ3=' + str(subst['delta3'])).center(len(headers[1]))
]
row += [(str(len(ds.valid_cases()))+'/'+str(ds.number_of_cases)).center(5+len(headers[2]))]
row += [str(len(number_oscillating[ds.name])).center(len(headers[3]))]
print ' '+' '.join([row[i] for i in range(len(row))])
Inspection of the table above shows 9 systems automatically identified with oscillatory potential, of which two are new design previously unidentified. In addition, 5 of these 9 designs have multiple phenotypes with oscillatory potential: S.4, S.10, S.11, S.12 and S.16. In the next subsections we focus our analysis on the new designs and explore the designs that have multiple potentially oscillatory phenotypes.
Analysis of newly identified designs with oscillatory potential
We analyze the phenotypes with oscillatory potential for the newly identified designs. Our aim is to determine if the phenotypes identified have dominant interactions that can be reduced to one of the original 7 designs identified in our previous analysis (10). This can be determined by inspecting the term signature or the S-System equations of the individual phenotypes with oscillatory potential.
We start with System S.3, by inspecting the S-System equations for the phenotype with oscillatory potential. To facilliate interpretation of the S-System equations, we remove the auxiliary variables by substituting their solution back into the dynamic equations.
ds = systems['S.3']
cases = ds([i for i in number_oscillating['S.3']])
for case in cases:
print 'Case ' + str(case.case_number) + ' [' + case.signature + ']'
ssystem = case.ssystem
# We remove the algebraic constraints by substituting their solution back into dynamic equations
ssystem = ssystem.remove_algebraic_constraints()
print 'Equations:'
display(ssystem.equations)
From these equations we can readily identify the effective modes of activator and repressor regulation. For Case 12, synthesis of \(x_1\) is regulated by a repressor-only mode of transcriptional control (only \(x_4\) influences synthesis); meanwhile, synthesis of \(x_3\) is regulated by a dual mode of transcriptional control (both \(x_2\) and \(x_4\) influence synthesis). Of the original 7 designs we analyzed, there are none that have this particular architecture; thus, this phenotype cannot be reduced to one of the original designs and hence this is a new design with potential to produce oscillations.
Next, we look at the S-System equations for the oscillatory phenotypes of design S.11,
ds = systems['S.11']
cases = ds([i for i in number_oscillating['S.11']])
print ds.name
for case in cases:
print 'Case ' + str(case.case_number) + ' [' + case.signature + ']'
ssystem = case.ssystem
ssystem = ssystem.remove_algebraic_constraints()
print 'Equations:'
display(ssystem.equations)
print ''
This design has two phenotypes with the potential to oscillate, Case 10 and Case 12. For Case 10, synthesis of \(x_1\) is regulated by a repressor-only mode of transcriptional control (only \(x_4\) influences synthesis); meanwhile, synthesis of \(x_3\) is regulated by an activator-only mode of transcriptional control (only \(x_2\) influences synthesis). This corresponds to the negative-only archictecture and thus can be reduced to the negative-only design [S.9]. For Case 12, synthesis of \(x_1\) is regulated by a repressor-only mode of transcriptional control (only \(x_4\) influences synthesis); meanwhile, synthesis of \(x_3\) is regulated by a dual mode of transcriptional control (both \(x_2\) and \(x_4\) influence synthesis). Again, this combination of dominant interactions cannot be reduced to those of one of the original designs.
The dominant interactions for Case 12 of S.3 and Case 12 of S.11 have the the characteristic architecture for their respective designs: repressor-only control of activator transcription and dual control of repressor transcription. Using this strategy, we have automatically identified both activator- and repressor-primary variants, where S.3 has a repressor-primary mode of repressor control (activator antagonizes repressor binding) and S.11 has a activator-primary mode of repressor control (repressor antagonizes activator binding).
Analysis of the design with four oscillatory phenotypes
Next, we focus our analysis on designs that have multiple oscillatory phenotypes to determine if these phenotypes can be co-localized in a 2-D slice of design space. We analyze the system with repressor-primary control of activator transcription and activator-primary control of repressor transcription (S.12) that has 4 distinct phenotypes with the potential to oscillate. We add constraints to the values of the \(\beta_2\) and \(\beta_4\) parameters for Case 18 to shift the automatically determined parameter values nearer to \(\beta_2\) = 1 and \(\beta_4\) = 1.
ds = systems['S.12']
psets=ds.co_localize_cases([i for i in number_oscillating['S.12']],['b2', 'b4'],
constraints=['$b2_1 > 100']
)
The values of \(\beta_2\) and \(\beta_4\) automatically determined for these co-localized phenotypes are shown in the table below.
print '\nTable of values for B2 and B4'
headers = ['Parameter'] + [('Case '+str(i)).center(12) for i in number_oscillating['S.12']]
print ' '+' '.join([i.center(9) for i in headers])
for name in ['b2', 'b4']:
row = [name]
for case_number in number_oscillating['S.12']:
row += ['%.3f'% psets[str(case_number)][name]]
print ' '+' '.join([row[i].center(len(headers[i])) for i in range(len(row))])
This slice is shown graphically by plotting the 2-D design space. We select the set of parameter values for the phenotype nearest to the middle (Case 18) and center the plot on the automatically determined parameter values for this case.
pvals = dspace.VariablePool(names=ds.independent_variables)
pvals.update(psets['18'])
fig = figure()
ax = pyplot.subplot(111)
cdict = ds.draw_2D_slice(ax, pvals, 'b2', 'b4',
[log(2), log(2)*1e6],
[log(2), log(2)*1e6],
intersections=[1, 3, 5])
The stability of the phenotypes is calculated, and shown graphically in the plot below. We overlay the boundaries of the different phenotypes and select four points within the regions of potential oscillation. The behavior of the full system will be simulated at the four points shown below.
fig = figure()
ax = pyplot.subplot(111)
points = [(5, 5, 'o'), (5, 3.5, 's'), (3.5, 0.6, 'v'), (0.8, -0.1, '^')]
pvals.update(psets['18'])
ds.draw_2D_positive_roots(ax, pvals, 'b2', 'b4',
[log(2), log(2)*1e6],
[log(2), log(2)*1e6],
color_dict=stability_colors,
resolution=300,
)
# We overlay the boundaries to help visualize the different phenotypes
ds.draw_2D_slice(ax, pvals, 'b2', 'b4',
[log(2), log(2)*1e6],
[log(2), log(2)*1e6],
colorbar=False,
facecolor='none', ec='k')
for i in points:
ax.plot(i[0], i[1], '.', marker=i[2], mfc='k', mec='k')
We test if the regions of potential oscillation result in sustained oscillations by simulating the full system. We define the original system, and process the strings so that they can be compiled by the python interpreter (we substitute '**' for '^').
original_system = ['x1. = B10*(a1/a10)*rho10^0.5*\
(rho1^-pi1 + (K20^g12)*delta1*x2^g12*K2^-g12\
+ (K40^g14)*rho1^-1*x4^g14*K4A^-g14)\
/(1 + (K20^g12)*delta1*x2^g12*K2^-g12\
+ (K40^g14)*x4^g14*K4A^-g14)\
- b1*x1 - mu*x1',
'xA. = BA0*(aA/aA0)*x1 - bA*xA - mu*xA',
'x2. = B20*(bA/bA0)*(n20*a10)*(K20^-1*B10^-1*rho10^-0.5)*xA - b2*x2 - mu*x2',
'x3. = B30*(a3/a30)*rho30^0.5*\
(rho3^-pi3 + (K20^g32)*x2^g32*K2R^-g32\
+ (K40^g34)*delta3*rho3^-1*x4^g34*K4^-g34)\
/(1 + (K20^g32)*x2^g32*K2R^-g32\
+ (K40^g34)*delta3*x4^g34*K4^-g34) - b3*x3 - mu*x3',
'xR. = BR0*(aR/aR0)*x3 - bR*xR - mu*xR',
'x4. = B40*(bR/bR0)*(n40*a30)*(K40^-1*B30^-1*rho30^-0.5)*xR - b4*x4 - mu*x4',
]
equations_rhs = [equation.split(' = ')[1] for equation in original_system]
equations_rhs = [equation.replace('^', '**') for equation in equations_rhs]
We define the dictionary of variable : dynamic equations that will be used for the CVODE interface.
equation_dict = {ds.dependent_variables[i]:equations_rhs[i] for i in range(len(equations_rhs))}
The original equations have been defined for the general class. We substitute the \(\delta_{\left\{1,3\right\}}\) and \(\pi_{\left\{1,3\right\}}\) parameters and the kinetic orders to obtain the equations for S.12.
subs = dict(substitutions[ds.name])
subs.update(g12=2, g14=2, g32=2, g34=2)
for key in equation_dict:
for old,new in subs.iteritems():
equation_dict[key] = equation_dict[key].replace(old, str(new))
Using these equations, we simulate the full system at the four points in the 2-D stability plot.
# initial_conditions = {i:100 for i in equation_dict}
pvals.update(psets['18'])
fig = figure()
fig.set_size_inches(6., 3)
times = [20, 20, 25, 20]
axes = [(-1.5, 0.5), (-1.5, 0.5), (-0.5, 1.5), (1,3)]
for i in range(len(points)):
pvals['b2'] = 10**points[i][0]
pvals['b4'] = 10**points[i][1]
params = dict(pvals)
params['mu'] = log(2) # Add mu from the original model.
case = ds(ds.valid_cases(p_bounds=pvals)[0])
initial_conditions = case.steady_state(pvals)
if i == 2:
initial_conditions['x4'] *= 10
# We remove the values for the auxiliary variables (not defined for the original systems).
_=initial_conditions.pop('X5')
_=initial_conditions.pop('X6')
t, x = integrate(equation_dict, initial_conditions, params, 0, times[i], 0)
ax=fig.add_axes([0.1+0.45*(i % 2), 0.1+0.45*(1 - i // 2), 0.3, 0.2])
ax.plot(t, np.log10(x['x4']), 'k')
ax.set_ylim(axes[i])
ax.set_yticks([axes[i][0]+0.5*j for j in range(int((axes[i][1]-axes[i][0])/0.5)+1)])
ax.set_xlabel('$t$')
ax.set_ylabel('$\log_{10}[x_4]$')
ax.plot(times[i]*0.05, (axes[i][1]-axes[i][0])*0.85+axes[i][0],
'.', marker=points[i][2],mfc='k', mec='k')
Simulation at the parameter values shows three phenotypes that produce sustained oscillations (Cases 16, 43 and 45) and one that produces damped oscillations (Case 18).
Conclusions
In this notebook we have shown the methods to analyze mathematical models of chemical and biochemical systems using the Design Space Toolbox V2. We have analyzed a total of 18 different models using an approach that does not require knowledge of parameter values, but uses automatically determined values and methods to explore the design space of these models. This approach enumerates the phenotypic repertoire of each of these systems to obtain a global perspective of the phenotypic potential. The behaviors we have identified are rich and include different signal amplification factors and dynamic characteristics. The methods used in this notebook are still under active development and are expected to be released in 2015 under the open-source GLPK V3 License.
References
Savageau MA, Coelho PMBM, Fasani RA, Tolla DA, and Salvador A (2009) Phenotypes and tolerances in the design space of biochemical systems. Proc. Natl. Acad. Sci. U.S.A. 106, 6435–6440.
Savageau MA, and Lomnitz JG (2013) Deconstructing Complex Nonlinear Models in System Design Space, in Discrete and Topological Models in Molecular Biology (Jonoska, N., and Saito, M., Eds.). Springer.
Lomnitz, JG, and Savageau, MA (2013) Phenotypic deconstruction of gene circuitry. Chaos 23, 025108.
Fasani RA, and Savageau MA (2010) Automated construction and analysis of the design space for biochemical systems. Bioinformatics 26:2601–2609.
Hunter JD (2007) Matplotlib: A 2D Graphics Environment, Computing in Science & Engineering 9, 90-95
Walt S van der, Colbert SC, and Varoquaux G (2011) The NumPy Array: A Structure for Efficient Numerical Computation. Computing in Science & Engineering 13, 22–30.
Pérez F, and Granger BE (2007) IPython: A System for Interactive Scientific Computing. Computing in Science & Engineering 9(3):21–29.
Oliphant TE (2007) Python for Scientific Computing. Computing in Science & Engineering 9(3):10–20.
Millman KJ, and Aivazis M (2011) Python for Scientists and Engineers. Computing in Science & Engineering 13(2):9–12.
Lomnitz JG, and Savageau MA (2014) Strategy Revealing Phenotypic Differences among Synthetic Oscillator Designs. ACS Synth Biol 3(9):686–701.
Vanderbei RJ (1996) Linear Programming: Foundations and Extensions (Springer, New York, NY).
Voit EO (2013) Biochemical Systems Theory: A Review. International Scholarly Research Notices 2013:e897658.
Atkinson MR, Savageau MA, Myers JT, and Ninfa AJ (2003) Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli. Cell 113(5):597–607.
Savageau MA (2010) Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology (CreateSpace Independent Publishing Platform). 40th Anniversary Edition.
Hindmarsh AC, and Brown PN (2005) SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software (TOMS) 31(3):363–396.