Background

Models of biological networks play important roles in our understanding of disease biology [22, 23], cancer [24], drug discovery [25], metabolic regulation [26], and many other subjects. However, simulation of large kinetic network models continues to be a major challenge, despite recent progress in high–performance simulation software [27–29]. The growth in size and complexity of biological pathway models has exceeded the growth of simulation hardware and software. In one study, a whole–cell M. genitalium model required 10 hours on a 128 node Linux cluster in order to simulate a single cell cycle [30]. Large–scale examples of kinetic simulations also arise in genome–scale kinetic models [31, 32]. Common simulation bottlenecks arise in parameter fitting and calibration of models, which require many simulations [33]. Thus, improvements in simulation performance are necessary for better and more comprehensive model fitting and to enable larger, more robust models.

In many real–world computing tasks, the relevant metric for performance is not the total computing power of the system, but rather the computations–per–watt. A compelling example is the adoption of specialized hardware designed for Bitcoin mining, which can easily exceed 40 times the performance–per–watt of a graphics processing unit (GPU) [34]. Cytomorphic chips operate at tens of milliwatts, yet in many cases still perform equal or faster simulations than desktop computers operating at tens of watts, representing a more than 1000–fold improvement in performance–per–watt for general networks [3, 6, 8, 9, 11]. This efficiency improvement may be used to package more units onto a die, thereby allowing more simulations to run in parallel. In addition, our current prototype hardware, based on a low–yield manufacturing process, is able to achieve up to a 30,000–fold speedup over Matlab and a 700–fold improvement over COPASI for stochastic reaction networks [6].

The present work focuses on generalized, digitally–programmable cytomorphic hardware described previously [3, 6, 8, 9, 11]. The hardware maps the thermodynamic laws that govern stochastic and deterministic molecular flux in chemical reactions to stochastic and deterministic electrical current flux in electronic transistors that also obey these same laws in a mathematically exact fashion. Thus, all biological model variables, including species concentrations, parameters, and reaction rates are represented in the hardware by electric current values. However, a naïve approach at solving ODEs in this way can lead to infinite integration errors due to parameter mismatches in analog circuits, also manifest in numerical integration errors on digital computers [6, 8, 9, 11] (it is hard to easily match analog parameters to more than 10 bits of precision on digitally programmable analog cytomorphic chips whereas digital components can routinely operate at 64–bit precision, even though they are modeling biological circuits that only operate at 2–5 bits of precision [4, 8, 9, 11]).

Key challenges

A major challenge in simulating biochemical models is ensuring flux symmetry (i.e., when species A is converted to species B, the rate of production of B should be exactly the rate of consumption of A). Consider the situation shown in Fig 1. This figure depicts a kinase cascade the active form of kinase A, represented as A^P, which in–turn phosphorylates kinase B, with their corresponding rates of change: (1)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Demonstration of the divergence problem.

In the phosphorylation cascade in (A), the quantities A + A^P and B + B^P should be constant in time. However, letting and results in the loss of this conservation relationship, as shown by the value of A^P in the numerical integration of this ODE system (B), which exceeds the total starting amount of A + A^P = 1. This phenomenon also applies to networks that do not have conserved quantities, as any steady–state value will tend to drift over time. Using the total quantity representation of Eq (2), this problem can be eliminated (C).

https://doi.org/10.1371/journal.pcbi.1008063.g001

Since the total amount of the kinase A is constant, the rate of change of A+ A^P should be zero. However, consider the case where perturbations are added to the kinetic constants for production of of A^P (letting and ) and B^P ( and ). It can be seen from these equations that this will only occur if either α and β are zero, or . However, the steady state values of A and A^P are already fixed by the equilibrium ratio of the first step in the cascade. Therefore, for general perturbations α and β the quantity A + A^P will change over time, violating conservation laws. In this example, we considered a kinase cascade because it shows how the divergence problem clearly violates conservation laws, but this phenomenon actually applies to all networks that reach a steady state, regardless of whether conserved quantities exist in the network or not.

To address this problem, the cytomorphic chip is designed to operate on conserved quantities of the system, which are mapped to conservation laws such as Kirchhoff’s current law in electronic circuits, a law that is always precisely obeyed [6, 8, 9, 11]. In the example in Fig 1, the conserved quantities are A + A^P = A_tot and B + B^P = B_tot. The system can be described using a pair of differential equations corresponding to A^P and B^P: (2)

Chip layout

Here, we briefly review the layout and specifications of the cytomorphic chip. A more complete description of the hardware can be found in [6, 8, 9, 11].

Fig 2 depicts the layout of the cytomorphic chip. The chip is composed of 20 blocks, each designed to solve a single biochemical reaction of the form: where the rate of this reaction is k_f * A * B − kr * C * D. In practice, modelers are accustomed to working with more complex reactions with lumped kinetic expressions such as Michaelis–Menten kinetics. However, physical processes at the molecular level invariably fall into this binary mass–action category. In protein complex formation, subunits are added one–at–a–time, and in enzyme catalysis, substrates bind the enzyme in an intermediate state. We will show later how lumped kinetic expressions can be used with the hardware.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Layout of the cytomorphic chip.

The chip is composed of 20 blocks of reaction units (A). Each block accepts as input total quantities for reactants A_tot and B_tot and outputs the total amount of product C_tot, as well as any remaining “free” reactant A_free and B_free which has not yet been converted into product. These current–based inputs and outputs are converted into digital signals before being sent outside the chip in order to prevent signal degradation.

https://doi.org/10.1371/journal.pcbi.1008063.g002

Terminology & validation methods

Table 1 lists definitions for terminology used in this article.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Terminology used in this article.

https://doi.org/10.1371/journal.pcbi.1008063.t001

In addition to producing programming files for the cytomorphic hardware, the cytomorphic compiler produces two other types of output that can be used to validate the compiler: (1) A Simulink model containing the blocks, parameters, and connections produced by the compiler, or (2) a differential equation system called a block simulation based on the block diagram of Fig 4. Either of these outputs can be used to simulate the circuit behavior over time, similar to the SPICE analog circuit simulator [38] used in circuit design, except that the simulations are transfer function–based (i.e. they are based on the block diagram of Fig 4, which uses gains, multipliers, and summation blocks for each stage instead of individual circuit components). In fact, these two formats are numerically equivalent, but serve different use cases. We use Simulink diagrams to visualize the block wiring, whereas we use block simulations to plot and compare compiler output in Jupyter notebooks. In the case of mass–action networks, the block simulation should correspond exactly to the SBML simulation. However, the underlying differential equations in the block simulation are based on total quantities, whereas SBML uses free quantities.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Simplified conceptual block diagram of a reaction unit (see [9, 11] for a detailed circuit–level diagram).

The chip has a total of 20 reaction units. Each reaction units block accepts as input total quantities for reactants A_tot and B_tot and outputs the total amount of product C_tot, as well as any remaining “free” reactant A_free and B_free which has not yet been converted into product. These current–based inputs and outputs are converted into digital signals before being sent outside the chip in order to prevent signal degradation.

https://doi.org/10.1371/journal.pcbi.1008063.g004

When evaluating the performance of the cytmorphic chip against conventional software simulation, it should be noted that some practical limitations exist in the current prototype of the hardware, which makes exact quantitative comparison with software simulations challenging. These include manufacturing variations in resistance / capacitance, the capacitance of the block integrator, analog–to–digital converter (ADC) clock speed and transistor mismatch. For variations that can be measured, we calibrate the compiler to adjust the parameters based on the magnitude of the variation. As described in detail later and as shown in S7 Fig, we can compensate for most of the output variation between the expected output and the output generated by the chip by adjusting the internal gains.

Despite the lower precision of analog simulations as compared to software simulation, it is often sufficient for running biological simulations. For example, Proctor et al. present simulations of the two stochastic models of the p53 signalling pathway [39]. The authors focus on whether their proposed models predict 1) the “existence” of sustained oscillations, 2) the shape of oscillations, 3) the effect of parameter changes or perturbations, and 4) the effect of stochasticity. In doing so, less attention is paid to the precise values of oscillation period, molecular concentrations, or parameters, partly because many of them are highly variable or not known. It is also because biological networks are often quite robust to variations in parameters.

Shift registers & block parameters

For the default configuration of the FF_EN_sw switches in Fig 4, it is apparent that the rate of production of the block’s main product, C, is given by: (3)

The forward and reverse rates for each block are determined by a combination of k_r and either KD_fw or KD_rv. Examining the equation shows that the forward and reverse rate constants are given by and respectively. Thus, the forward and reverse dissociation constants can be used to tune the relative forward and reverse rates, whereas the overall reaction rate of all blocks in the network can be tuned by changing the k_r value. Also, any block that produces a species can also serve as a degradation reaction by utilizing ratC.

SRAM & network building blocks

The chip’s SRAM is used to program connections between the input and output ports of the blocks. While each protein block is designed to simulate a fundamental second-order biochemical reaction, we can capture more complex dynamics by connecting several blocks together. Since it is instructive to visualize how total variables affect block connection in non-trivial ways, we demonstrate topologies involving multiple blocks that simulate a 3-stage feed-forward network cascade and a “fan-out” reaction. These reactions are chosen because they form the basis for most larger networks and exhibit several recurrent wiring topology that appear in more complex building blocks. These building blocks are selected from elementary types of subnetworks—they are not related to motifs (which are defined by statistical over–representation in real networks [40, 41]). Instead, building blocks are comprised of simple subnetworks such as “fan–in” and “fan–out” configurations, where a single species is produced or consumed by multiple reactions respectively. For a more in-depth discussion on these circuit networks, we refer interested readers to previous publications [8, 9].

Fig 5 is a graphical abstraction of a single protein block. The currents flowing into the input ports represent the concentration of reactants in a forward reaction. The currents flowing out of the output ports represent the concentration of products formed from the reaction. Since each reactant and product may participate in several downstream reactions, the concentration currents for Afree, Bfree and Ctot are copied using current mirrors. In what follows, we use species names S, T, etc. to distinguish species in the network from the port names on the block (A, B, C, and D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Input and output ports for the cytomorphic block.

The two main inputs A_tot and B_tot are used to compute the forward rate. Internally, the block subtracts its own C_tot value from A_tot and B_tot to compute A_free and B_free (see Fig 2B). The chip’s main output is C_tot. If there is another reaction that consumes C, then C_tot should be connected to the A_tot or B_tot input of the consumer block. Otherwise, C_tot should be connected to the block’s own C_free input to allow calculating the reverse rate. If the reaction is reversible, D_free may also receive input from another block or be simply wired to unity in the case of a single product C. Other ports describing the blocks forward and reverse rates are used in in certain building blocks described below. The chip has 2 copies of Afree, Bfree, fw_up, rv_up, and 5 copies of Ctot (2 negative and 3 positive). In addition, the chip also copies the values of its C_free and D_free inputs to the Cfree_cp and Dfree_cp output ports respectively. These copies are used to route the input values to additional blocks in fan–in configurations (see below).

https://doi.org/10.1371/journal.pcbi.1008063.g005

To simulate a 3-stage feed–forward biochemical cascade, we connect three protein blocks to simulate a single irreversible production reaction, two reversible substitution reactions, and an irreversible degradation steps as shown in Fig 6. We need just three protein blocks instead of four to simulate these reactions because each block is designed to simulate both a second order reaction and product degradation. In this case, block 3 is wired to simulate both the substitution reaction between T and U and the degradation of U. To represent the production of S and its conversion to T, we connect Ctot of block 1 to Atot of block 2. This connection is repeated between block 2 and 3 for the substitution reaction of T to U. To simulate the reverse reaction fluxes, we connect output ports that represent the concentration of free products to input ports Cfree to produce the right concentration of free species. Finally, to account for the decrease in total variable of species S and T as U degrades, we backpropagate the degradation flux of U upstream by connect ports rv_up from block 3 and 2 to ports Cdeg from block 2 and 1 respectively. These design principles hold true for larger feed-forward circuit cascade as well.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. A three–step feed–forward network.

This network is comprised of a linear chain of unimolecular mass–action processes (A). To convert this network to a wiring, we designate blocks 1–3 as the main produces of species S, T, and U respectively. Block 1 produces S, hence its main output port C_tot is connected to the A_tot input of block 2 and so on. This signal is summed with the initial value of T (if non–zero). Since the cytomorphic chip uses currents for computation, summing signals is achieved simply by connecting multiple signals to the same input port. The last block in the chain, which produces U, also serves as a degradation reaction . The block has its ratC parameter set to the degradation rate k_d. Additionally, the amount of U degraded must also be subtracted from the total values of S and T. The rv_up port computes the total loss in C for each block and is propagated to the block immediately upstream by connecting to the Cdeg port. The Antimony/SBML model for this network (B) is converted to a block simulation that gives identical output to the SBML simulation (C). Since this is a mass–action network, numerical differences can be made arbitrarily small by adjusting integrator tolerances. All connections created by the compiler are shown in the wiring diagram. More info is available at https://github.com/cytocomp1/cytocomp/tree/master/case-studies.

https://doi.org/10.1371/journal.pcbi.1008063.g006

In the “fan–out” reaction shown in S1 Fig, two protein blocks were used to model the two reversible reactions. Similar to the feed-forward cascade, ports Cfree on block 1 and 2 are wired to ports Ctot so that the free concentration of U and W can be used to simulate the correct reverse flux. Since there is no degradation reaction, we do not need to backpropagate any loss of product via rv_up. Instead, to calculate the unbound concentration of S, we simply subtract U and W from input which is achieved by closing an internal use-it-and-lose-it feedback loop in block 1 and adding a negative current Ctot of block 2 respectively. We have included two additional examples for a dissociation reaction and a “fan–in” reaction in the supporting information (S2 and S3 Figs).

Lumped kinetics

The term “lumped kinetics,” as used here, refers to any kind of process in a model that represents multiple elementary steps, where an elementary step is defined as a bimolecular mass action reaction (e.g. substrate binding / unbinding, or the catalytic step in enzyme catalysis). A common instance of this is enzyme kinetics: E + S ⇌ ES → E + P, which represents an enzyme E that converts substrate S into product P. These reactions occur at the following rates: (4)

The individual forward and reverse rate constants of the binding step are difficult to measure directly, so these two processes are usually condensed into a single process by assuming either rapid equilibrium of the binding process (which leads to the well–known Michaelis–Menten kinetics [42]) or by assuming the enzyme–substrate complex is at steady–state (Briggs–Haldane kinetics [43]). In either case, the rate law for the resulting lumped process can be expressed as: where for Briggs–Haldane kinetics or k_r/k_f for Michaelis–Menten kinetics.

Our general approach to simulating these lumped expressions on the cytomorphic chip is to break them down into the constituent steps of Eq 4. However, this re-creates the problem of determining the forward and reverse rate constants k_f and k_r respectively. From the lumped expression, we can determine k_cat and K_M. However, we need one additional constraint to specify both k_f and k_r. This constraint comes from the upper limit of the chip’s simulation speed.

Consider the block diagram of a reaction unit containing a negative feedback loop around the part of the circuit that processes A as highlighted in Fig 7. In designing electronic amplifiers, it is common to account for the so–called phase margin. In an amplifier circuit as well as in the highlighted feedback loop, there exists the possibility that the output signal can be 180° out–of–phase with the input. Since the feedback is negative, the signal will be inverted and cause constructive interference with the input, leading to instability. Furthermore, this feedback loop also possesses a parasitic pole due to the current mirror that produces the A_free signal. Taken together, these conditions lead to the stability rule [11]: where C_par is the parasitic capacitance at the gate node of the A_free current mirror, C = 0.1μF is the capacitance of the integrator capacitor, and r is the overall rate of the block (used as a scaling factor for both the forward and reverse rates). In the preceding expression, we assumed that A_free < B_tot. The roles of A and B can be reversed if this is not the case. This equation can be simplified to: (5)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Negative feedback around the circuit for A.

https://doi.org/10.1371/journal.pcbi.1008063.g007

This gives an upper bound for the value of k_f based on the global value and the local value which depends on the reaction and simulation conditions. These two values can be condensed into a single constant ρ, called as the margin, which in general varies per reaction. We will encounter several more examples of lumped expressions that expand into elementary processes that occur at rate ρ. In practice, we can assume a reasonable lower bound for ρ and refine the estimate after simulating the network. If we assume C_par = 1pF, this gives an upper–bound of . The value of is model and simulation–dependent, but by assuming a reasonable upper bound of 100, we obtain ρ = 1000. This value can be used to run digital simulations for a given hardware configuration which in turn can be used to refine the value ρ.

Returning to enzyme kinetics, the higher the value we choose for k_f (and hence k_r), the more rapid the enzyme–substrate binding. Since Michaelis–Menten kinetics are derived based on an equilibrium assumption, a larger k_f will tend to make this assumption more valid. Therefore, choosing yields the best approximation to Michaelis–Menten kinetics without causing instability.

Fig 8 shows the result of plotting the dynamics of the elementary binding / catalysis network versus the original Michaelis–Menten lumped single–process network for various values of ρ.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Different margin values and their respective simulations.

In these simulations, E = 1, S_initial = 10, k_cat = 1, and K_M = 10. For the elementary two–process network, the reverse binding rate is calculated automatically by the compiler using the supplied value of ρ as k_r = ρK_M − k_cat if or k_r = ρK_M otherwise.

https://doi.org/10.1371/journal.pcbi.1008063.g008

It is worth observing that, from a design standpoint, our compiler uses the two–step network to “approximate” models with lumped processes such as Michaelis–Menten kinetics. However, from a mechanistic standpoint, Michaelis–Menten kinetics represents an approximation of the corresponding physical two–step process. In effect, this line of reasoning has taken us full-circle from a mechanistic representation to a lumped process and back again. However, for the purpose of practical modeling, this is a necessary detour, since it is feasible to measure the K_M and V_max values for enzymes, but in general it is not possible to measure binding apart from catalysis (i.e. the parameters k_f, k_r, and k_cat).

The meaning of this observation is that while elementary, mass–action processes are biologically accurate, they are not feasible from a modeling standpoint. Thus, modelers should continue to use kinetics that can be parameterized with quantifiable parameters, while specialized hardware should continue to use elementary processes that are conducive to efficient implementation. Breaking–down these high–level expressions into low–level expressions, what we hereafter refer to as “expansion”, is one of the main functions of the cytomorphic compiler.

A limitation of this “expansion” method is that when the substrate is not in excess of the enzyme (i.e. E ≪ S does not hold), then there is significant deviation between the two–step mechanistic process and the idealized Michaelis–Menten approximation (the mechanistic process will tend to lag behind the lumped process), regardless of the value of ρ. However, Michaelis–Menten kinetics also relies on the assumption that E ≪ S, and thus would be a physically inaccurate modeling assumption in this case.

Matching algorithm for lumped kinetics

In order to successfully expand lumped kinetic expressions into constituent components, it is necessary to (1) identify, from the rate law, what type of lumped expression is represented, and (2), obtain the values of all lumped constants and use these to compute the individual rate constants.

One approach to solving (1) would be to simpify and canonicalize the rate law expression and compare this simplified version with each known kinetic expression on a tree–basis (we refer to canonical expressions like the Michaelis–Menten formula as “archetypes”). However, this approach is sensitive to different factorizations of the expression and requires that all archetypes also be in canonical form. Instead, we have opted to use an algorithm for determining the equivalence of expressions based on hash–coding [44].

In 1971, William Martin observed that, by assigning numeric values to symbols appearing in a mathematical expression, an equivalence class can be computed for a given expression by simply evaluating it on a finite field. Two mathematical expressions can thus be evaluated for equivalence by simply comparing their “hash codes,” so–constructed. Based on this observation, we have constructed an algorithm that operates on a database of common rate–laws and their “hash codes.” We call these common rate–law expressions archetypes, and pre–compute their hash codes according to the algorithm below. We then compute these pre–computed values to input expressions to determine which type of lumped kinetics a reaction employs. The algorithm is split into two stages (lines 1–12 and 14–25 in Fig 9 respectively).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Pseudocode for the Martin matching algorithm.

https://doi.org/10.1371/journal.pcbi.1008063.g009

The first stage accepts as input an SBML kinetic law expression in the form of an abstract syntax tree (AST). The leaves of this tree may be either the substrate / product species that participate in the reaction, or symbols that refer to numerical constants (such as K_M). The algorithm assigns a different prime number to each symbol in the AST if it is a substrate / product of the reaction (Fig 9 line 7) and unity otherwise. Symbols that are assigned relatively prime values may be used in finite field expressions (such as addition, multiplication, etc.) and the results will still contain relatively prime factors, which helps to distinguish the hash code values after many operations.

On the other hand, when constructing archetypes, any symbolic constant (such as K_M) can in principal contain any value and any symbolic name in a given SBML model. Hence, we assign a value of 1 to these constants. Hence, we must also assign a value of 1 to all constants in the matching phase.

The first stage then computes the hash code of the expression to the known values for all archetypes in the database according to Martin’s algorithm (Fig 9 line 9). The output of the first stage is the matched archetype along with the positions of the substrates and products in the expression (Fig 9 line 12).

Once the input expression has been mapped to a valid archetype, we can proceed to map kinetic constants within the expression (for example, to sort out which symbols in the input expression correspond to K_M versus V_max). For each constant, we assign a numeric value, setting all other constants to unity (Fig 9). The assigned constant, as well as all substrates and products, must be relatively prime. We then iterate through all SBML parameters (Fig 9 line 21) and substitute the same prime number for the value of this parameter in the input expression 9 line 23). If the computed hash of the input expression matches the archetype, the parameter p in the input expression corresponds to the symbol s in the archetype and we identify the role of p as being either K_M, V_max, etc. in the returned datastructure.

It is evident from Fig 9 that the complexity of the first stage of the matching algorithm is O(N_sN_a), where N_s is the number of symbols in the input expression, and N_a is the number of archetypes in the database. Fortunately, most kinetic models make use of a limited repertoire of kinetic laws, such as Michaelis–Menten, Botts–Morales, and Monod–Wyman–Changeaux kinetics. This limits both N_s and N_a in most practical models.

The second stage of the algorithm has O(N_pN_a) complexity, where N_p is the number of symbols in the input expression which map to SBML parameters, and N_a is the number of symbols in the archetype a. For the same reasons as above, these are also capped in most practical models.

Gene regulatory kinetics

Another major type of lumped kinetics occurs in models of gene regulatory networks. Consider the LacI repressor, which controls expression of the lac operon (lacO) in bacteria. The LacI repressor is a homotetramer, but might be better described as a dimer of dimers. Each dimer subunit contains a DNA–binding site for lacO. The repressor binds to lacO as a two–step process. Binding of allolactose to LacI causes the repressor to enter an inactive state P (protruded) with decreased overall affinity, releasing the operator site. Using lumped kinetics, the transcription rate of the operator is given by [45, 46]: where ν is the transcription rate, α is an experimentally determined rate constant, K_m is the equivalent equilibrium constant of the two binding steps, R is the concentration of the active repressor, and α₀ is the basal level of transcription under fully repressed conditions (the “leakage” rate).

In order to expand this process to a simulatable form, we could decompose binding into a two–step process. However, the cytomorphic chip provides a Hill function that allows this two–step process to be modeled as a single step (labeled n in Fig 2). Using the Hill function, we could write the overall binding process as: where R is the active repressor, O is the operator site, and B is the bound (repressed) complex. Transcription can then be modeled as a simple first–order process without substrate depletion (both the A_FB_EN and B_FB_EN in Table 2 should be off). This expanded model exhibits a time delay proportional to τ = kr + (k_f ⋅ R)² compared to the lumped expression, which assumes rapid equilibrium of the binding process. Our approach to minimizing this discrepancy is the same as in the enzyme kinetics case—we maximize the forward rate constant k_f subject to the margin ρ defined in Eq 5. The expanded model can reproduce the dynamics of complex networks to a high degree of accuracy, as shown below.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Tunable parameters for each block.

https://doi.org/10.1371/journal.pcbi.1008063.t002

In a highly influential study in 2000, Elowitz et al showed that a genetic oscillator (the “repressilator”) can be constructed from three genetic repressors [45]. To validate the cytomorphic compiler’s ability to translate repressor kinetics, we obtained a dynamical model of the repressilator from the BioModels database [47, 48]. This model contains a total of 12 reactions, half of which are degradation reactions. To reduce the number of blocks required to encode the model, we condensed the degradation reactions into the production rate laws for the three genes and proteins in the system: (6) (7) and similarly for Y and Z. We constructed a new SBML model based on these simplifications and generated a block configuration as shown in S4 Fig. The configuration is largely analogous to the three-stage feed-forward cascade presented earlier and repeated three times for the three gene and protein species. They differ mainly in how they account for the total variable. In the case of the repressilator, transcription factors and mRNAs regulate transcription and translation but they are not directly converted to their downstream molecules the way S is converted to T and U in the feed-forward cascade. The total variable for each biomolecule in the repressilator is not dependent on their downstream products and consequently, the degradation fluxes are not backpropagated up the circuit, as seen in the unconnected rv_up ports.

Fig 10 shows a digital SBML simulation of the repressilator model, a block simulation (emulation of the cytomorphic hardware using a digital computer), and the cytomorphic chip data for this model. The compiler is successful in configuring the wiring between blocks to produce the oscillations seen in the repressilator. The difference in oscillation amplitudes from the hardware simulation can be attributed to manufacturing variations in the chip. We hypothesized that these variations could be explained in terms of variation of the internal gains within the chip. The internal gains are the parameters ratC, KDfw, KDrv, kr, and kdeg in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Comparison of repressilator model simulations.

Also shown is an SBML simulation of the original model (A), a block simulation (B), and data collected from the cytomorphic chip (C). Due to manufacturing variations, blocks in the cytomorphic chip have different gains, hence the peak heights are different. The time axis on the chip data plot corresponds to “simulation time,” i.e. the actual duration of running the simulation, as opposed to “model time”, the time according to the dynamics of the model. To be useful, a hardware simulation should take less time to run than the timescale according to the model’s dynamics. This allows for multiple in silico expriments to be performed in for the amount of time a single physical experiment would take. To quantify this ratio, we performed correlation analysis on the chip data against the block simulation based on the wiring on S4 Fig (D). The cross correlation shows a peak at ≈ 6.5 seconds, indicating a six–fold difference between model and simulation time. This does not represent a speedup over a software SBML simulation, but the cytomorphic chip exhibits constant scaling up to the maximum number of blocks.

https://doi.org/10.1371/journal.pcbi.1008063.g010

In order to compensate for these variations, we sought to determine if the gains in the digital block simulation could be adjusted to match the chip data. To perform these adjustments, we used the Nelder–Mead method with an objective function of the mean–squared–error (MSE) between the block simulation and the chip data for mRNA and protein levels. In practice, we found it sufficient to adjust only the KDfw and KDrv gains, which enabled us to reduce the MSE from 124 nA to 15.6 nA, an 87% decrease. The results of this adjustment are plotted in S7 Fig. Indeed, this strategy has been previously exploited in analog hardware for cochlear processors to correct for variation and enable deaf subjects to successfully use cochlear prosthetics on the first try [3, 49, 50].

Higher–order compilation

While the requirement that all processes in a model be reducible to mass–action kinetics may seen restrictive, there is a very large class of models that consist only of this type of process. The field of rule–based modeling is a very active area of research (see [51] for a review). Rule–based models are composed of multi–state species. For example, a protein can have multiple phosphorylation sites, and each site can exist in either a “phosphorylated” or “unphosphorylated” state. These rules can be used to generate a network of all enumerable molecular states, or alternately simulated stochastically without enumeration using agent–like methods such those employed by the simulator NFSim [52].

Rule–based models invariably generate mass–action networks when enumerated or otherwise reduced to a simulatable form (this is not strictly a requirement, but non–mass action networks are a rare use–case for rule–based modeling and we do not consider them here). Furthermore, whereas NFSim scales linearly with the number of rules [52], the cytomorphic chip has constant scaling up to the maximum number of blocks (with the ability to connect to additional chips and thus increase the maximum network size in the future).

Taken together, these factors suggest that the ideal application of cytomorphic hardware could be the simulation of rule–based mass–action networks. We therefore selected a model of a MAPK signaling cascade from the rule–based modeling platform PySB [53] and based on a previous study [54] to validate our cytomorphic compiler. This cascade consists of the MAPK ERK and its upstream activators Ras and Raf. A contact map for this model, generated using RuleBender [55], is shown in Fig 11.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. A rule–based MAP kinase model [53, 54] and its corresponding block configuration produced by the cytomorphic compiler.

(A) A contact map for the kinase cascade generated using RuleBender [55]. The contact map shows the molecular species present in the model—Ras, Raf, MAPK/ERK Kinase (MEK), extracellular signal-regulated kinase (ERK), MAP kinase phosphatase (MKP), and Protein phosphatase 2 (PP2A). When expanded to an SBML mass–action network representation, this rule–based model expands into a network with 20 reactions and 21 distinct dynamical states. An SBML simulation of the flattened rule–based model (B) and a block simulation (C) are identical up to 3 significant figures (simulated using libRoadRunner, CVODES solver, absolute and relative tolerances 10⁻²⁰ and 10⁻¹² respectively).

https://doi.org/10.1371/journal.pcbi.1008063.g011

The full block layout and wiring for this model is shown in S6 Fig.

Rule–based models are one approach to managing complexity. They allow the user to specify models in terms of biomolecules with multiple states (such as multiple phosphorylation sites and multiple binding sites that can either be occupied or not) and automatically enumerate all possible discrete states. Similarly, electronics designers have long used automated logic synthesis to generate chip layouts from high–level logic specifications. Just as automated placement and routing was a necessary technology for enabling rapid growth in complexity of integrated circuits (ICs), we believe technologies such as the cytomorphic compiler presented here will be necessary for rapid growth of biomimetic electronics.

Discussion

Special–purpose hardware is necessarily designed to solve a specific subset of problems. One must sacrifice generality for the sake of improved performance and / or efficiency on this subset. Our system clearly makes this trade–off of generality for performance. In this section we discuss limitations, work–arounds, and future plans for improving the generality of the hardware and software.

The most obvious limitation of the cytomorphic hardware is that all processes are constrained to bimolecular mass–action reactions. However, a large body of rule–based models exists which are not hindered by this restriction, since rule–based models can be converted into mass–action reaction networks [51]. We have also provided an automatic method for reducing enzyme and repressor binding kinetics to mass–action networks. Nevertheless, our method requires these lumped kinetic expressions to be explicitly built–in to the compiler. Our future plans include adding support for user–defined lumped kinetic expression reductions to allow for kinetic expressions not already handled.

Other less–obvious limitations include dynamic range and forcing function support. The cytomorphic hardware exhibits less dynamic range than digital simulations because the state variables in the cytomorphic chip are physical currents. Very large values would tend to be attenuated due to limitations in the amount of current that can be used by the chip. In most cases, the overall concentration values used in the model can be rescaled without changing the model dynamics by also rescaling kinetic constants appropriately. As for forcing functions, SBML allows model quantities to be defined in terms of user–specified, potentially time–dependent expressions (called assignment rules). These rules can be used to implement time–varying input to the system. In control theory, such inputs are called forcing functions. This feature in SBML allows researchers to test their models against different types of time–varying input. However, the cytomorphic hardware itself cannot synthesize arbitrary waveforms, thus preventing this feature from being used. Instead, assignment rules are evaluated at the initial state of the model. Any subsequent changes are not accounted for. In theory, this could be remedied by connecting an arbitrary waveform generator to the cytomorphic hardware. We plan on exploring this possibility as the hardware matures. A similar argument applies to other SBML input–centric features (e.g. events and rate rules).

One application where a hardware implementation of an ODE solver is particularly useful is in estimating the parameters of a biological model from known experimental data. Current methods of parameter estimation iterate between model simulation using an initial set of parameters, and exploring a new set of parameters that minimizes a loss function [56]. It is widely recognized that the bottleneck in parameter estimation is the numerical integration of the system of ODE since it is repeated for every new set of parameters. Furthermore, a robust statistical approach to parameter estimation must necessarily make use of multiple solutions of an optimization problem, either as part of a Bayesian approach (e.g. [57]) or uncertainty analysis [58]. Thus, a highly parallel, high performance–per–Watt analog computer for performing many such simulations is highly sought–after.

In contrast to digital simulations, which are exact, analog simulations suffer from accuracy problems that are caused by mismatches in component resistances and capacitances due to manufacturing variations. We have shown that these variations can in principle be compensated for by adjusting the internal gains of each reaction block in the chip (S7 Fig). We hope to further develop this procedure to engineer a robust platform for parameter fitting based on analog computation.

The cytomorphic chip is analog hardware that runs in continuous-time domain. It contains noise generators that amplify thermal noise to create thermodynamically-accurate random fluctuations for biochemical reactions. These factors enable the chip to simulate multiple stochastic reactions in parallel to enjoy network–size–invariant speedups. By running the simulation on the cytomorphic chip, we can leverage the chip’s speedup to reduce the time of each solution to the parameter fitting problem [8, 9]. In more concrete terms, our current prototype consists of 20 reaction units and can simulation roughly 2 ~1–second repressilator simulations simultaneously. In contrast, the digital compute we used consumes roughly 100 Watts of power and can run approximately 30 analogous simulations per second. The cytomorphic chip can be scaled up to ~3000–fold more reaction units before reaching the same thermal envelope as the workstation, which is roughly a 200–fold improvement in performance for the same thermal envelope. Attaining this limit would require addressing significant challenges in data input / output and signal fidelity, among other considerations, but suggests that analog hardware may become an attractive platform for simulation in the future.

Disclaimer

The content is solely the responsibility of the authors and does not necessarily represent the views of the National Institutes of Health, the Betty Moore Foundation, or the Alfred P. Sloan Foundation.