Computational design of molecular motors as nanocircuits in Leishmaniasis

Two strategies for the computational design of the nanocircuit

A coupled bistable switch and repressilator, with the bistable switch under the control of repressilator can be delivered using a liposome (of nano size) with antibodies (Ab) targeting the lipophosphoglycan (LPG) displayed on an L. major infected macrophage. The liposome will be designed to withstand the acidic conditions of the phagolysosome of targeted macrophages. The circuit with bistable switch under the control of repressilator is given in Figure 1.

Figure 1. Nanocircuit encapsulated in an Ab coated nanoliposome.

The repressilator has TetR (represses P_LTet01) under the control of the P_Llac01 promoter, lacI (represses P_Llac01) under the control of λP_R, λcI (represses λP_R), which is under the control of P_LTet01. All the repressilator proteins repress each other in a cyclic fashion. The bistable switch comprises of cIts (represses P_Ls1con), which is under the control of P_Llac01 and lacI (represses P_Llac01) under the control of P_Ls1con. P_Llac01 is repressed by lacI of both the repressilator and bistable switch. The reporter gene is placed downstream of P_Llac01.

The addition of an inducer of PLTet01, namely aTc, will increase the production of λcI, which will inhibit λP_R. As production of lacI of the repressilator is stopped, some of the repression on P_Llac01 will be released. cIts represses the production of lacI of the bistable switch. Thus, the reporter gene present downstream of P_Llac01 will be expressed in increased levels.

The expression will be kept under control as the repression of lacI production in repressilator will release the repression on TetR (represses P_LTet01), which will counteract with the inducer aTc.

The choice of the reporter gene in the nanocircuit may be based on one of the following strategies.

Strategy I (Figure 2). The reporter gene may be devised such that its transcription results in ssRNA (antisense RNA), which will be complementary to the mRNA strand coding for the amino acids 809–823 – LQWVEEASNMFPDF (Batters et al., 2014). The inhibition of this calmodulin binding (CB) motif in Myosin XXI will result in the hindrance in CB, necessary for the motility of the parasite. There will be an increase in the levels of calmodulin calcium complex in the system. The calcium calmodulin complex activates the myosin heavy chain kinases, which will result in the continuous contraction of myosin. Thus, the motility and intracellular trafficking of the parasite will be disturbed. Also, the complex of calcium and calmodulin will choose to bind to other CB proteins, namely calcineurin. Activation of calcineurin will in turn activate PP1. PP1 and PP2A function as molecular switch that reciprocally regulate the eukaryotic phosphatases, hence suggesting the fact that PP1 modulates PP2A activity in cells. (Ohki et al., 2001). As more PP1 is activated, phosphorylation of ABC transporter proteins necessary for its regulation will also be affected, thus disturbing the efflux of molecules. Consequently, L. major will be put under stress as the motility, efflux of endogenous metabolites; survival mechanisms and intracellular trafficking are affected, leading to death of the parasite.

Figure 2. Flow chart of strategy I.

CBF, Calmodulin Binding Motif; PP1, protein phosphatase 1; ABC transporters, ATP-binding cassette transporters.

Strategy II (Figure 3). The reporter gene can be devised to synthesize Protein Phosphatase 2A (PP2A), which dephosphorylates the myosin heavy chain. Excess production of PP2A by the circuit will result in continuous relaxation of myosin. This will in turn affect the parasite’s motility and intracellular trafficking, endangering its survival. PP2A will also disturb the phosphorylation states needed for ABC transporter activity (Berridge MJ, 2010). Also, dephosphorylation of mitogen-activated protein kinase (MAPK), by PP2A, has been noted in several studies (Toivola & Eriksson, 1999). Several studies have reported that since extracellular signaling kinase (ERK) is involved in activation of transcription factors, like Ap-1, which is associated with ABCC1 (multidrug resistance protein-1) expression, dephosphorylation by PP2A may also have an effect on ABC transporter expression (El Azreq et al., 2012; Xu et al., 2006). Thus, the efflux of molecules will be affected. The trafficking mechanism of L. major, its survival mechanisms and drug resistance property breaks down as a result. Thus, L. major will cease to exist, unable to cope up with the build-up of stress.

Figure 3. Flow chart of strategy II.

MAPK, Mitogen activated protein kinase; ABCC1, ATP Binding Cassette Subfamily C Member 1; PP2A, Protein phosphatase-2A; Ap-1, Activator protein 1; ABC transporters, ATP-binding cassette transporters

Cofactor – streptavidin and biotin. Streptavidin has been selected as a cofactor to ensure the internalization of the reporter produced by the circuit. This is based on the fact that biotin-streptavidin-based isolation of Leishmania has been widely practiced. Biotin, if encapsulated inside the nanoliposome containing the nanocircuit, will become attached to the surface proteins of L. major inside the macrophages. Therefore, the streptavidin-reporter complexes attaches with this bound biotin, finally leading to internalization of complex (Figure 4).

Figure 4. Use of Biotin-streptavidin for internalization of reporter.

Designing and simulation of the nanocircuit using Tinker Cell

The nanocircuit was designed using different biological parts available in Tinker Cell (http://www.tinkercell.com/). The transcriptional repression processes were also indicated in the circuit. Sequences for the different parts of the gene modules were obtained from various sources: the reporter DNA and protein sequences were obtained from NCBI; the DNA and protein sequences for repressor proteins, ribosome binding sites, terminators were obtained from Registry for Standard Biological Parts; and the promoter sequences with binding strengths for the repressor proteins were acquired using TRANSFAC and TFBind (Supplementary File 1). These sequences were entered in appropriate places in the Text Attributes dialog box in Tinker Cell. Steady state parameter scan was performed for the designed circuit. Simulation of the working of the circuit was done by altering the values of protein degradation rates, rbs and promoter strengths, dissociation constants, and hill coefficients, until the desired graph was obtained. Inducer was added to PLTet01 of the circuit, and the simulation was performed again. The final results were exported in SBML format for further analysis. The SBML files were validated using an online SBML validator (http://sbml.org/Facilities/Validator/).

Validation of the nanocircuit using Genetic Regulatory Network Inference using Time Series (GRENITS) and BoolNet

The SBML files of the simulated circuit (with and without inducer) were loaded in COPASI (v4.18; software for simulation and modeling of biochemical networks; http://copasi.org/) and time course data for a period of 10s was generated for the protein concentrations. The acquired files were analyzed by GRENITS package (v1.24.0; https://www.bioconductor.org/packages/release/bioc/html/GRENITS.html) for the circuit’s convergence and network link probabilities. Inferred network of the circuit components was obtained from the network links with probability greater than 0.8. The files were also analyzed by BoolNet (v2.1.3; https://cran.r-project.org/web/packages/BoolNet/index.html). The possible attractor states formed by the circuit, the probability of Boolean network transitions and network wiring were found. Robustness of the circuit for around perturbations was also checked. Attractor results from BoolNet were exported as .net format and the circular layout of the state transitions were obtained using Pajek (v4.10; http://mrvar.fdv.uni-lj.si/pajek/).

Studies on circuit behavior using Berkeley Madonna

The SBML files of the designed and simulated circuit (with and without inducer) were loaded in COPASI, and the ODE files of the circuit were exported in .mmd format. The ODE files were loaded to Berkeley Madonna (v8.3.18; www.berkeleymadonna.com). The Bistable switch - Repressilator coupling equation was included as follows:

                cod1 = LcI1+ (PLlac01_strength*rs4_c*rs5_c * DefaultCompartment)

The ODEs were integrated using Euler’s method from 0 to 10s. A plot of Antisense RNA or PP2A vs time was obtained. The plots for the circuit with and without inducer were compared to determine the type of gene expression in the nanocircuit and to understand the effect of inducer on the response time delay. A plot of Antisense RNA or PP2A vs LacI3, showing the phase plane of the two proteins, was obtained.

A nullcline was plotted using the option Nc in Berkeley Madonna. Another nullcline was obtained by plotting several trajectories from different initial conditions using the option Ic.

The point of intersection of the two nullclines shows the steady state of the designed circuit.

Design and simulation of the nanocircuit

The circuit was designed using Tinker Cell, a CAD software. All the sequences for the different parts of the gene modules were entered as detailed in the Methods section.

Figure 5 shows the designed nanocircuit without inducer.

Figure 5. Design of nanocircuit.

The upper portion of the circuit forms the repressilator with TetR, Lambda cI, LacI protein repressing each other in a cyclic fashion. The lower portion forms the Bistable switch comprised by cIts and lacI mutually repressing each other. LacI of repressilator represses PLLac01 of Bistable switch forming the coupling between repressilator and Bistable switch. The reporters Antisense RNA or PP2A, as well as streptavidin, are placed downstream of P_LLac01 of Bistable switch.

Simulation was performed and the desired graph was obtained for the values shown in Table 1.

Figure 6 shows the working of the circuit in the absence of the inducer aTc. cIts and so Antisense RNA or PP2A (and streptavidin) are in off stage, due to the repression by LacI2, while LacI3 is on. The repressilator proteins TetR, LcI, LacI2 are forming oscillations.

Figure 6. Simulation of working of the nanocircuit.

Figure 7 shows the circuit with inducer aTc added, which will reduce the repression caused by TetR on P_LTet01.

Figure 7. Design of circuit with inducer.

Table 2 shows the values of the three parameters of the induced activation, for which the desired graph was obtained. The rest of the values of the circuit were kept constant. The graph (Figure 7) shows the working of the circuit after the addition of aTc.

The reporters (Antisense RNA or PP2A) are switched on when Antisense RNA or PP2A degradation rate is in the range of 0–2s, while the lacI3 gets switched off, as in Figure 8.

Figure 8. Simulation of working of the nanocircuit after addition of inducer.

The cofactor Streptavidin also peaks in the same range (Figure 9).

Figure 9. Behavior of Streptavidin after addition of inducer.

The repressilator proteins TetR, LcI and LacI2 form limit cycle oscillations in this range (Figure 10).

Figure 10. Oscillation of repressilator proteins after addition of inducer.

Thus, the circuit was designed and the working in the absence and presence of inducer. The corresponding values of the parameters were also noted.

Validation of the nanocircuit

The validation of the circuit with and without inducer was performed using GRENITS and BoolNet package in R (v3.2.2; https://cran.r-project.org/).

Convergence of the circuit was checked using Monte Carlo Markov Chain Algorithm by GRENITS. Various parameters, such as indicator variables of Gibbs sampler, coefficient of regression, network connectivity parameter, precision of each regression and intercepts of regression of Chain 1 and Chain 2 of state transitions were analyzed.

The graphs in Figure 11 shows the convergence of indicator variables of Gibbs sampler (Gamma), coefficients of regression (B), network connectivity parameter (Rho), precision of each regression (Lambda) and intercepts of each regression (Mu) of state transition chain 1 and 2, thus validating the circuit without inducer.

Figure 11. Convergence plots of the nanocircuit without inducer.

Graphs in Figure 12 shows that convergence of the parameters was also obtained for the circuit with inducer, thus validating it.

Figure 12. Convergence plots of the nanocircuit after addition of inducer.

The links between the different genetic components in the circuit were analyzed using GRENITS package. Links with probability > 0.8 were considered to infer the network.

Figure 13 shows the link probabilities between the components. From the plot, it is evident that LacI3 of the bistable switch is the most important and strong regulator of the circuit, as it regulates the expression of all other proteins.

Figure 13. Link probabilities of the nanocircuit.

Figure 14 shows the number of regulating parents for each of the genes and the probabilities of each of the parent regulating it. It shows that LacI3 has no parent regulating it.

Figure 14. Network uncertainty of the nanocircuit.

The inferred network between the genetic components considering the links with probabilities > 0.8 (Figure 15A) is shown in Figure 15B

Figure 15.

(A) Link probability threshold and (B) the inferred network for the nanocircuit.

Similarly, links in the nanocircuit with inducer were analyzed using GRENITS, which resulted in the same results; LacI3 was the main regulator (Figure 16) and had no parents regulating it (Figure 17).

Figure 16. Link probabilities of the nanocircuit with inducer.

Figure 17. Network uncertainty of the nanocircuit with inducer.

The inferred network from the links with probability > 0.8 (Figure 18A) is shown in Figure 18B.

Figure 18.

(A) Link probability threshold and (B) the inferred network for the nanocircuit with inducer.

The inferred network shows that all the components are linked to one another with high probability. This proves the non-randomness of the network.

The possible attractor states formed by the circuit were determined using BoolNet to visualize the switching on and off of the genes comprising the circuit.

Figure 19 shows the ten possible attractor states with continuous switching of the repressilator genes between on and off states. The reporters are off when LacI3 is on and they are on when LacI3 is off, which is an evidence for the toggling between the genes of the bistable switch.

Figure 19. Attractor states of the nanocircuit.

Figure 20 depicts the ten attractors with either 8 or 16 states forming the basin of attractor.

Figure 20. State graph of the nanocircuit.

The attractor state analysis of the nanocircuit with inducer also yielded 10 attractors with 8 or 16 states forming the attractor basin (Figure 21 and Figure 22).

Figure 21. Attractor states of the nanocircuit with inducer.

Figure 22. State graph of the nanocircuit with inducer.

Figure 23 illustrates the Probability Boolean Network transitions in the circuit with and without inducer. It gives the binary representation of the attractors and its basins.

Figure 23. Probability Boolean Transitions.

The robustness of the circuit to 10 different perturbations was checked using BoolNet (a state transition table was generated, one or several transitions are perturbed randomly and the gene transition function are robust from the modified transition table).

Figure 24 shows the percentage of attractors present in the perturbed networks of the circuit: 10% of original attractors were present in 6 perturbed networks; 20% in 1; 30% in 2; and 40% in 1 network. Hence, the result showed that 70% of the random results had the original attractors undisturbed.

Figure 24. Percentage of original attractors in perturbed networks (nanocircuit).

Figure 25 shows the Gini index of state-in degrees of the circuit. The state-in degrees of a particular state is the number of transitions leading to it. The Gini index of the in-degrees is returned as a characteristic value of the network. It is a measure of inequality. If all states have an in-degree of 1, the Gini index is 0. If all state transitions lead to one single state, the Gini index is 1 (Kim et al., 2007). The higher the Gini index, the higher the number of transitions leading to it. This parameter is measured based on the fact that in biological networks, many state transitions lead to the same states. Figure 25 show that 96% of the states have a high Gini index, proving that the designed circuit is a biologically valid network and not a random network.

Figure 25. Gini index of state-in degrees (nanocircuit).

Similar results were obtained for nanocircuit with inducer with BoolNet. In total, 70% of the perturbed networks had the original attractors (Figure 26), while 96% of the attractors had high Gini index of state-in degrees (Figure 27). Hence, the biological validity of the circuit has been proved.

Figure 26. Percentage of original attractors in perturbed networks (nanocircuit with inducer).

Figure 27. Gini index of state-in degrees (nanocircuit with inducer).

Based on the transition function of each gene found from the links between the gentic components, network wiring was made using BoolNet. From the network (Figure 28), it is evident that LacI3 controls the expression of other proteins of the circuit and is also a function of itself, while not being controlled by other components.

Figure 28. Network wiring in the nanocircuit.

Figure 29A and B shows the circuit layout of the 2¹²⁹ transitions and the probable attractors in the circuit and the circuit with inducer, respectively.

Figure 29.

Circular layout of the state transitions of the nanocircuit (A) without and (B) with inducer.

Thus, several other relationships are possible that can lead to a difference in network behaviour. Here the concept of modularity is of paramount importance. Modularity is the degree to which a system’s components can be segregated or recombined. In synthetic biology, each module has a fixed arrangement pertaining to which it exhibits a characteristic behaviour; if this modular arrangement is altered there are perturbations in the global dynamics, oscillations and kinetics of a system. A second round of simulation performed without parameter scan sketched 16 attractor states which portrayed a stochastic behaviour of attractor probabilities. The state matrix which we framed from the graph of transition states made it is pretty much clear that both proteins decay over the same time period in our analysis when the circuit is initially inactive (no stimuli) and then switch to the active state. This behavioural features work in a contrast fashion when the system switches in an opposite direction i.e., from initially active to inactive state. Our Boolean analysis solidified the idea that below a threshold value, expression is never induced. The next set of quantitative analysis which was performed using the GRENITS module gave results which were congruent with the qualitative justifications.

Studies on circuit behavior

The plot of Antisense RNA or PP2A vs time (Figure 30 and Figure 31) was obtained using Berkeley Madonna. The gradual increase in the concentration of the reporter shows that there is graded expression in the circuit. Comparison between the two graphs shows that in addition to a high level production of reporter, the response time delay has also been reduced by the addition of inducer.

Figure 30. Response time delay in nanocircuit.

Figure 31. Response time delay in nanocircuit with inducer.

Determination of steady state

Considering the biological fact that at equilibrium protein concentration does not change over time, nullclines (dx/dt = 0 or dy/dt = 0) have been plotted for reporters and LacI3 in the phase plane of the two proteins, i.e. d(Antisense RNA or PP2A)/dt = 0 or d(LacI3)/dt=0). The point of intersection indicates the steady states (d(Antisense RNA or PP2A)/dt = 0 and d(LacI3)/dt=0)). The nullclines of a bistable switch genes intersect thrice with two stable steady states. From Figure 32 and Figure 33, it is evident that the designed nanocircuit favors only one steady state with high reporter concentration and low LacI3 concentration. The individual trajectories in the graph show the evolution of the system from their respective initial conditions.

Figure 32. Geometric nullclines of the nanocircuit.

Figure 33. Geometric nullclines of the circuit with inducer.

From the sharp threshold of steady state in equilibrium conditions, it is apparent that the circuit is ultrasensitive.

To evaluate the dynamics and evaluation of the circuit under non-equilibrium conditions, quasipotential landscape was mapped. It also gives a quantitative measure of the epigenetic landscape. The directionality of the circuit can also be studied.

Figure 34 shows the contour plot of quasipotential energy at different levels of reporter and LacI3. The regions marked by yellow lines are local minima in the quasipotential surface. It corresponds to the stable states as the areas of minimum potential that have the highest probability of occupancy. This plot also shows that the circuit works by autoactivation. Autoactivation of the reporter leads to amplification of transient differences in expression between the reporter and LacI3.

Figure 34.

Contour plot of the Quasipotential of nanocircuit (A) without and (B) with inducer.

The 3D quasipotential landscape is given in Figure 35. The branching valleys and ridges depict the stable cellular states and the barriers between those states, respectively. The “deeper” valleys in the figure are associated with a higher probability of occupancy than the “shallower” valleys. Figure 35 shows that autoactivation has caused multistability in the circuit. A circuit with an initial condition of low reporter and high LacI3 will flow down-hill to the intermediate stable state with minimal quasi potential (Figure 36). For the circuit to flow to the stable steady state from the intermediate state with high reporter and low LacI3, it has to cross the potential barrier between the two states, for which additional reactions are required which will need energy. Therefore, the probability of this transition is less. In Figure 36, it can be seen that there is a smooth down-hill flow to the stable steady state, with high reporter levels and low LacI3 levels from all initial conditions when inducer is added. There is translational burst of the reporter, leading it to the dominant attractor. Therefore, the probability of occupancy of this stable steady state is very high. The number of intermediate stable states with minimal potential has been reduced by the inducer.

Figure 35.

3D Quasipotential landscape of the nanocircuit (A) without and (B) with inducer.

Figure 36.

Linearized quasipotential landscape of the nanocircuit (A) without and (B) with inducer.

Figure 36 shows the linearized 3D plot from which the relative distances between the stable states can be inferred.

In the absence of inducer, there is a higher number of intermediate attractors closely spaced (Figure 36A). When inducer is added, the dominant attractor is placed well away from the initial condition with a smooth down-hill flow (Figure 36B). This occurs due to the high level expression of the reporter placed downstream of the active promoter and low levels of regulatory LacI3, which is upstream of this promoter.

Thus, it can be concluded that the circuit works by autoactivation of the reporter, while repressing LacI3. This has led to multistability. Addition of inducer reduces the multistability caused by the autoactivation in the circuit and smooths the undulating quasipotential landscape. Therefore, the circuit will occupy the stable state of high reporter levels with high probability.