pystablemotifs - Control tutorial#
import pystablemotifs as sm
import networkx as nx
from timeit import default_timer
Model importation#
We show how to use biolqm to import a model and convert it to BooleanNet format:
import biolqm
lqm = biolqm.load("http://ginsim.org/sites/default/files/2018_zanudo_proliferation.zginml")
biolqm.save(lqm, "Proliferation.txt", "booleannet")
'Proliferation.txt'
Load network and find attractors#
See the Basic Usage Tutorial for further details.
primes = sm.format.import_primes('Proliferation.txt',remove_constants=True)
sm.format.pretty_print_prime_rules(primes)
AKT* = PIP3
E2F* = !Rb
EIF4F* = mTORC1
FOXO3* = !MAPK | !AKT
GFs* = GFs
MAPK* = RAS | PIP3
PI3K* = RTK | RAS
PIP3* = PI3K
PRAS40* = !AKT
Proliferation_b1* = Proliferation_b1 & Proliferation_b2 | EIF4F & S6K | E2F
Proliferation_b2* = E2F & EIF4F & Proliferation_b1 & S6K
RAS* = RTK
RTK* = GFs & !MAPK & !S6K | FOXO3 & GFs
Rb* = !cycE
S6K* = mTORC1
TSC* = !AKT & !MAPK
cycE* = AKT & !FOXO3 | E2F
mTORC1* = !PRAS40 & !TSC
ar = sm.AttractorRepertoire.from_primes(primes)
ar.summary()
There are 3 attractors.
{'AKT': 0, 'E2F': 1, 'EIF4F': 0, 'FOXO3': 1, 'GFs': 0, 'MAPK': 0, 'PI3K': 0, 'PIP3': 0, 'PRAS40': 1, 'Proliferation_b1': 1, 'Proliferation_b2': 0, 'RAS': 0, 'RTK': 0, 'Rb': 0, 'S6K': 0, 'TSC': 1, 'cycE': 1, 'mTORC1': 0}
{'AKT': 0, 'E2F': 0, 'EIF4F': 0, 'FOXO3': 1, 'GFs': 0, 'MAPK': 0, 'PI3K': 0, 'PIP3': 0, 'PRAS40': 1, 'Proliferation_b1': 0, 'Proliferation_b2': 0, 'RAS': 0, 'RTK': 0, 'Rb': 1, 'S6K': 0, 'TSC': 1, 'cycE': 0, 'mTORC1': 0}
{'AKT': 'X', 'E2F': 1, 'EIF4F': 'X', 'FOXO3': 'X', 'GFs': 1, 'MAPK': 'X', 'PI3K': 'X', 'PIP3': 'X', 'PRAS40': 'X', 'Proliferation_b1': 1, 'Proliferation_b2': 'X', 'RAS': 'X', 'RTK': 'X', 'Rb': 0, 'S6K': 'X', 'TSC': 'X', 'cycE': 1, 'mTORC1': 'X'}
Define a control target#
Select a set of node values that we wish to drive the system toward. In this example, we specify a set of nodes (of size 1), namely Proliferation_b1=1
, that identifies two attractors. The succession-based methods require that the target is consistent with at least one attractor.
target = {'Proliferation_b1':1}
Search for knockins/knockouts that achieve the target#
Brute-force#
The max_drivers
parameter limits our search to a maximum number of concurrent interventions.
Note that the brute force approach scales poorly with the size of the network and unless a value is specified for every variable, it does not guarantee convergence to an attractor. Therefore, the intervention must be permanent in general.
start=default_timer()
interventions = sm.drivers.knock_to_partial_state(target,primes,max_drivers=2)
end=default_timer()
print("Time running method:",end-start)
print("Sets found:")
for x in interventions:
print({k:v for k,v in sorted(x.items())})
Time running method: 0.03742701800001669
Sets found:
{'AKT': 1}
{'PI3K': 1}
{'Rb': 0}
{'mTORC1': 1}
{'RTK': 1}
{'RAS': 1}
{'cycE': 1}
{'PIP3': 1}
{'E2F': 1}
{'Proliferation_b1': 1}
{'MAPK': 1, 'PRAS40': 0}
{'GFs': 1, 'MAPK': 0}
{'PRAS40': 0, 'TSC': 0}
{'FOXO3': 1, 'GFs': 1}
{'EIF4F': 1, 'S6K': 1}
Grasp search#
Here we use a heuristic approach to search for drivers that achieve the target. The GRASP_iterations
parameter controls how many heuristic searches are performed.
GRASP_iterations=2000
start=default_timer()
interventions = sm.drivers.GRASP(target,ar.primes,GRASP_iterations)
end=default_timer()
print("Time running method:",end-start)
print("Sets found:")
for x in interventions:
print({k:v for k,v in sorted(x.items())})
Time running method: 0.15597427600005176
Sets found:
{'PIP3': 1}
{'cycE': 1}
{'RAS': 1}
{'PI3K': 1}
{'RTK': 1}
{'AKT': 1}
{'E2F': 1}
{'Rb': 0}
{'mTORC1': 1}
Internal history#
In this method, all succession diagram pathways that are consistent with the target are identified. At each branch point in the succession diagram, the desired target stable motif is searched for internal driver node sets that drive the system into a narrower trap space containing the target. All possible paths are considered. All interventions can be permanent or temporary. Convergence to a consistent attractor (if it exists) is guaranteed.
start=default_timer()
interventions = ar.reprogram_to_trap_spaces(target,
target_method='history',
driver_method='internal')
end=default_timer()
print()
print("Time running method:",end-start)
print("Sets found:")
for x in interventions: print({k:v for k,v in sorted(x.items())})
Time running method: 0.00047217600058502285
Sets found:
{'GFs': 1}
{'GFs': 0, 'Rb': 0}
{'E2F': 1, 'GFs': 0}
{'GFs': 0, 'cycE': 1}
Minimal history#
This method also selects drivers for target stable motifs at each succession diagram branch point. It differs from the previous method in that it does not require these drivers to all be internal to each stable motif. This allows the method to uncover more parsimonious interventions at the cost of increase computational burden. It may identify interventions that are inconsistent with the target; such interventions must be temporary (e.g., temporary administration of a drug).
start=default_timer()
interventions = ar.reprogram_to_trap_spaces(target,
target_method='history',
driver_method='minimal')
end=default_timer()
print()
print("Time running method:",end-start)
print("Sets found:")
for x in interventions:
print("---")
print("One temporary intervention from each list, in order.")
print("("+str(len(x))+" interventions in total)")
for y in x: print(y,"\n")
Time running method: 0.0025028489999385783
Sets found:
---
One temporary intervention from each list, in order.
(1 interventions in total)
[{'GFs': 1}]
---
One temporary intervention from each list, in order.
(2 interventions in total)
[{'GFs': 0}]
[{'Rb': 0}, {'E2F': 1}, {'cycE': 1}]
GRASP history#
This method is like the two above, but the driver search is conducted using a heuristic approach. This is most useful in extremely large networks. The benefit of the GRASP method is that it does not consider all possible variable combinations, and can therefore consider larger driver sets with comparitively little additional computational burden.
start=default_timer()
interventions = ar.reprogram_to_trap_spaces(target,
target_method='history',
driver_method='GRASP',
GRASP_iterations=500)
end=default_timer()
print()
print("Time running method:",end-start)
print("Sets found:")
for x in interventions:
print("---")
print("One temporary intervention from each list, in order.")
print("("+str(len(x))+" interventions in total)")
for y in x: print(y,"\n")
Time running method: 0.15195031399980508
Sets found:
---
One temporary intervention from each list, in order.
(1 interventions in total)
[{'GFs': 1}]
---
One temporary intervention from each list, in order.
(2 interventions in total)
[{'GFs': 0}]
[{'cycE': 1}, {'Rb': 0}, {'E2F': 1}]
Minimal merge#
In this method, all minimal trap spaces containing only attractors consistent with the target are found, and a brute force search is conducted to identify interventions of minimal size that drive the system into these trap spaces. Unlike the brute-force method, it does not require that the intervention be permanent. Interventions that are inconsistent with the target must be temporary. Generally, this method is slower than others, but also finds the most parsimonious interventions. The worst-case computation time grows rapidly with the max_drivers
parameter, as all possible sets of variables up to this size can be considered.
start=default_timer()
interventions = ar.reprogram_to_trap_spaces(target,
target_method='merge',
driver_method='minimal',
max_drivers=4)
end=default_timer()
print()
print("Time running method:",end-start)
print("Sets found:")
for x in interventions: print(x)
Time running method: 0.021936680000180786
Sets found:
{'GFs': 1}
{'GFs': 0, 'PI3K': 1}
{'Rb': 0, 'GFs': 0}
{'GFs': 0, 'RTK': 1}
{'GFs': 0, 'RAS': 1}
{'GFs': 0, 'cycE': 1}
{'GFs': 0, 'PIP3': 1}
{'GFs': 0, 'E2F': 1}
Internal merge#
This method is like the above, but it only considers interventions that are internal to the stable motifs that constitute the trap spaces under consideration. Typically, this is faster, but it has the potential to overlook some interventions. Interventions can be temporary or permanent. As with the previous method, the worst-case computation time grows rapidly with the max_drivers
parameter; however rather than considering combinations of all variables, this method considers only combinations of variables that belong to the stable motifs that make up the target trap space. Therefore, the scaling is better than the previous method.
start=default_timer()
interventions = ar.reprogram_to_trap_spaces(target,
target_method='merge',
driver_method='internal',
max_drivers=4)
end=default_timer()
print()
print("Time running method:",end-start)
print("Sets found:")
for x in interventions: print(x)
Time running method: 0.000921599999855971
Sets found:
{'GFs': 1}
{'cycE': 1, 'GFs': 0}
{'Rb': 0, 'GFs': 0}
{'E2F': 1, 'GFs': 0}
GRASP merge#
This method is like the two above, but the driver search is conducted using a heuristic approach. This is most useful in extremely large networks when it is anticipated that only large intervention sets will drive the system to its desired target. This is due to the fact that method’s computational cost scales polynomially with the size of considered intervention set (whereas the minimal merge method scales combinatorially).
start=default_timer()
interventions = ar.reprogram_to_trap_spaces(target,
target_method='merge',
driver_method='GRASP',
GRASP_iterations=500)
end=default_timer()
print()
print("Time running method:",end-start)
print("Sets found:")
for x in interventions: print(x)
Time running method: 0.13632410999980493
Sets found:
{'GFs': 1}
{'GFs': 0, 'PIP3': 1}
{'GFs': 0, 'PI3K': 1}
{'GFs': 0, 'RAS': 1}
{'GFs': 0, 'Rb': 0}
{'cycE': 1, 'GFs': 0}
{'GFs': 0, 'RTK': 1}
{'E2F': 1, 'GFs': 0}