8  Gene regulation in time

8.1 Introduction

In this practical you will get hands on experience into making gene regulatory networks models. First you will practice how to encode regulatory interactions into ordinary differential equations and logical rules. Next you will study how many interactions produce self-sustained activity configurations (attractors) in a Boolean network model (Rozum et al. (2023)). Lastly, you will learn how to predict the genes that cause cell fate transitions (attractor changes).

Today we will use the following python file: circuitsfunctions.py. In it you will find the code of the functions we use throughout the practical, have a quick look at it.

VS Code

If you have been experiencing problems with Spyder and restarting the kernel, we suggest you install VS Code - this is an integrated development environment where the code runs faster, and kernel issues should be gone. Follow this installation instructions.

8.2 Regulatory interactions

Open Boolean_practical.py and circuitsfunctions.py. The first file has the main running code, and the second one has only the functions that are used in the main.This is a way to keep the code “simple” and easy to read. First have a look at the functions ODEgeneRegulation() and logicalRule(); look at what each takes as input and returns as output.

Both functions model an \(AND\) gate where the variable \(A\) and variable \(B\) positively regulate \(C\) This could represent that \(A\) and \(B\) are transcription factors that form a protein complex, and that it is via this complex that they can bind the promoter of \(C\) to promote its expression.

While ODEgeneRegulation() encodes this regulation with an ordinary differential equation that requires several parameters, logicalRule() only needs the logical operator relating the input genes.

ODE implementation

\(\frac{dC}{dt}=p\frac{A^n}{K^n+A^n}.\frac{B^n}{K^n+B^n}-dC\)

Where \(p\) is the production rate, \(K\) is the half saturation constant (represents the affinity of a TF to the promoter), \(n\) is the cooperativity of transcription factor binding, and is the same for \(A\) and \(B\), \(d\) is the decay rate.

Logical rule implementation

\(C_{t+1}=A_t\) & \(B_t\)

Python code for ODE and logical rule

def ODEgeneRegulation(a,t,parameters): 
    prod=parameters['prod']
    decay=parameters['decay']
    Ksat=parameters['Ksat']
    variableA=parameters['variableA']
    variableB=parameters['variableB']
    n=parameters['n']
    outputC=a[0]
    doutputC=prod*variableA**n/(Ksat**n+variableA**n)*variableB**n/(Ksat**n+variableB**n)-decay*outputC #Tip: Change this in exercise 8.3
    return(doutputC)

def logicalRule(variableA,variableB):
    return(variableA and variableB)

Exercise 8.1 (Mathematical thinking) What is the minimum information you need to encode a regulatory interaction with either function? In what cases would you prefer to use an ODE or a logical rule for a model?

Let’s run the model to simulate what happens if \(A\) and \(B\) are both active/expressed. Do this using the ODEgeneRegulation() and logicalRule() functions. Copy the following lines of code inside the main in Boolean_practical.

• Let op! Running the code may show a NumbaDeprecationWarning, ignore it.

# ODErun has three arguments: model, geneA, geneB
A=10
B=10
ODErun(ODEgeneRegulation,A,B) 

# Boolean model
# look at the text in the terminal for the result:
A=1
B=1
print("the boolean operation of variableA ",A," AND variableB",B," is:", logicalRule(A,B))

Both implementations of this regulatory interaction produce the same effect, i.e. \(C\) is active when both \(A\) and \(B\) are active. Yet, the ODE version recovers quantitative differences in \(C\) activity with values higher than 1.

Now, let’s explore the output of each function using different values of \(A\) and \(B\). We are going to use the code below to plot the output of the ODE and the Boolean model in a 2D heatmap. Notice that the exploration values for the ODE model is 11 values (0,1,2,3,4,5,6,7,8,9,10), while for a Boolean model is 2 (0 or 1). Copy this code inside the main in Boolean_practical.py.

ODE model simulation - AND gate:

explorationvalue = 11 # an ODE model allows us to explore more values than a Boolean model
ode_output = np.zeros((explorationvalue, explorationvalue))
for variableA in range(0, explorationvalue):
    for variableB in range(0, explorationvalue):
        parameters = {'prod': 0.01, 'decay': 0.001,'Ksat': 4, 'n': 2,'variableA':variableA,'variableB':variableB} # prod, decay, Ksat, n, and initial values for A and B
        cells = odeint(ODEgeneRegulation, 0, np.arange(0, 1000.1 , 0.1), args=(parameters,)) 
        ode_output[variableA, variableB] = cells[-1, 0]

Boolean network simulation - AND gate:

explorationvalue=2 # a Boolean model assumes there is only 2 possible states: 0 or 1
bool_output = np.zeros((explorationvalue, explorationvalue))
for variableA in range(0, explorationvalue):
    for variableB in range(0, explorationvalue):
        bool_output[variableA, variableB] = variableA and variableB #AND #Tip: Change this in exercise 8.3

After running each, this function plots the ODE and Boolean results side by side:

ODEBooleanPlot(ode_output, bool_output)

Change the running time for the ODE model, and the \(K\) of each input, and see what effect that has on \(C\).

Exercise 8.2 (Biology) What similarties and differences in the activity of \(C\) do you see using each approach? Why is the activity of \(C\) higher in the ODE output than the Boolean one?

Exercise 8.3 (Altorithmic thinking:) Now let’s compare how ODE and Boolean models represent different regulatory interactions. You can use the same code in the previous box. Remember to modify the logical operator in the Boolean section, and to modify the \(doutputC\) equation in the ODEgeneRegulation function in circuitsfunctions.py. As before change the parameters in the ODE to see how that changes the activity of \(C\):

1. Either A or B can activate C.
2. A represses C.
3. A and B together repress C.
4. A represses C and B activates it.
5. A or B can activate C, but not when both are active at the same time.
6. Include an extra input variable, and make a multi-dimensional logical gate, or try any other biological scenario yon can think of.

Exercise 8.4 (Biology and Algorithmic thinking) Is there a logical gate that can be better represented with an ODE than with Boolean logic? Give an
examples of biological regulatory interactions that can be described with each logical gate?

8.3 Gene regulatory network

Now let’s move to a network model made of many individual regulatory interactions (Garcı́a-Gómez et al. (2020)). This model includes many regulatory interactions discovered for the cells of plant roots through years of experimental research. Here, all this knowledge is encoded as logical rules.

The model consists of 18 variables representing transcription factors, hormones, peptides. Look at the rootNetwork() function in the circuitFunctions.py file. This function defines the state of 18 variables (CK, ARR1, SHY2, AUXIAAR, etc.) based on an input stored in parameters, then it updates the state of each variable using the logical operators \(AND\), \(OR\) and \(NOT\); the result is stored in w_variable. The function returns the update state w_variable for each of the 18 variables.

Time course simulation

Let us define a random initial condition for each of these 18 variables, and the total number of time steps we will solve the network using the logical functions. First, we need to input the initial condition of the variables to rootNetwork() using parameters. Next, we save the output of this function as the new current state. We repeat this procedure iteratively for as many time steps as we defined, saving in a matrix the time course of the simulation. This can be done using a for loop as shown below:

timesteps=20
variables=18
matrix = np.zeros((timesteps+1, variables), dtype=int) 
matrix[0,:] = np.random.randint(0, 2, size=variables) #Random initial condition
for i in range(timesteps):
    # we pass the current state
    parameters= {'CK': matrix[i,0], 'ARR1': matrix[i,1], 'SHY2': matrix[i,2], 'AUXIAAR': matrix[i,3], 'ARFR': matrix[i,4], 'ARF10': matrix[i,5], 'ARF5': matrix[i,6], 'XAL1': matrix[i,7], 'PLT': matrix[i,8], 'AUX': matrix[i,9], 'SCR': matrix[i,10], 'SHR': matrix[i,11], 'MIR165': matrix[i,12], 'PHB': matrix[i,13], 'JKD': matrix[i,14], 'MGP': matrix[i,15], 'WOX5': matrix[i,16], 'CLE40': matrix[i,17]}
    # we save the new state
    matrix[i+1, :] = rootNetwork(parameters) # we save the output in i+1

plotBooleanTimecourse(matrix,timesteps)

At the end the function plotBooleanTimecourse() will make a plot to see how each variable in the network changes ts activity throughout the simulated timesteps. 100 timesteps is enough to reach an attractor for this network.

• For the initial condition you can also use any combination of 0 and 1 you like, or one of the attractors reported in the paper and check what happens when you update them.

Exercise 8.5 (Algorithmic thinking) Try this initial condition:

matrix[0,:]=[0,1,1,0,1,1,0,0,1,0,1,1,1,1,1,0,1,0]

Solve this initial condition using rootNetwork() and rootNetworkAsynchronous(). Compare the output in each case. What happens and why? What does an asynchronous update do?

This network has 18 variables, and then \(2^{18}\) = 262,144 possible states. We can solve each of these conditions (and wait…), or instead explore just a subset of them…

Exercise 8.6 (Algorithmic thinking) Using the previous code, add for loop to solve 100 random initial conditions, and then save the final activity configurations (attractors) in an attractors matrix.

The attractors matrix should have the dimensions (initial_conditions, variables) so that you can plot the results using the plotBooleanAttractors() function.

• Notice here you should save only the final state for each initial condition you explore.

• Tip: 100 initial conditions is a good number to explore.

ICs=100
attractors = np.zeros((ICs, len(parameters))) 

# your code

plotBooleanAttractors(attractors) # it takes as argument a matrix of attractors

You should have found many different attractors. Let’s use a multidimensional reduction technique to see how these attractors relate to each other. For this you need to install the umap-learn package using this line:

pip install umap-learn

More details on how to install packages here: Introduction to python.

Now that you have installed umap-learn, let’s run a UMAP of the attractors you recovered (this takes a bit of time to run, so be patient 😉):

UMAPBoolean(attractors)

Exercise 8.7 (Biology / Conceptual thinking) Q5 Why do you see clearly defined groups and not a continuous distribution of attractors?

Now let’s analyse the attractors you found. First, let’s group similar attractors using the sorted() function:

attractors_sorted = np.array(sorted(attractors.tolist()))
plotBooleanAttractors(attractors_sorted) 

Exercise 8.8 (Biology / Conceptual thinking) Q6 Do some attractors occur more frequently than others. Why is this?

Exercise 8.9 (Biology) Some variables are active (1) in the majority of the attractors, others in half, and others in very few. In how many attractors is the variable SHR active? What does this suggest about its regulation?

Now let’s remove the repeated rows (duplicate attractor states) to see unique attractors using the np.unique() function:

_, unique_indices = np.unique(attractors, axis=0, return_index=True)
attractors_unique = attractors[np.sort(unique_indices)]
plotBooleanAttractors(attractors_unique) 

Exercise 8.10 (Biology / conceptual thinking) How many unique attractors did you find? Compare them with the ones reported in the paper. Are they all fixed points?

8.4 Cell differentiation – jumping from one attractor to another

Depending on what we want to answer a continuous or discrete model may be more appropriate for a model. To study the role of many genes in cell differentiation, a Boolean model might be better particularly if we lack details of the parameters underlying each reaction. Then if we want to study how cells jump from one state to another, a continuous model might be more appropriate.

Here we will see how we can convert the Boolean model to a continuous one and use it to predict which regulators are able of causing a change in the state of the system (changes in cell fate!)

Compare the code of the functions rootNetwork() and rootNetworkODE(). Notice how in rootNetworkODE() the logical rules are represented with \(min\) and \(max\) functions (explained in the lecture), and then used in a sigmoidal function to create a continuous ODE model.

• \(AND\) operator is a \(min\) function, \(OR\) operator is a \(max\) function, and \(NOT\) operator is \(1-x\).

Use the code below to run a random initial condition for the 18 variables (IC) and see how the system behaves.

timerunning=10.1 
times = np.arange(0, timerunning, 0.1)

IC = np.random.randint(0, 2, size=18).tolist() #random initial condition
parameters = {'decayrate': 1, 'h': 50} 
cells = odeint(rootNetworkODE, IC, times, args=(parameters,)) 

plotODEroot(cells,times)

Exercise 8.11 (Modeling choices) Do the attractors match those recovered with the Boolean network?

Now let’s use the model to study cell differentiation routes. Let’s start in the following initial condition (IC vector), and find a change in a variable that produce a jump to another attractor (end). For this you can simply flip the activity of a gene in the initial condition (from \(0→1\) or \(1→0\)) and see the final attractor matches the desired end state. This simulated change mimics how an environmental signal may impinge in the regulatory network.

# We start here: 
IC=[0,0,0,0,1,0,1,1,1,1,1,1,1,0,1,0,1,0]
# We want to end here. 
end=[1,1,0,0,1,1,1,1,1,1,0,0,1,0,0,0,0,1]
# variable order
# CK, ARR1, SHY2, AUXIAAR, ARFR, ARF10, ARF5, XAL1, PLT, AUX, SCR, SHR, MIR165, PHB, JKD, MGP, WOX5, CLE40

# your code 

plotODErootTransition(cells,times) # use this function to plot your result (it already knows the desired end state)

Exercise 8.12 (Algorithmic thinking / biology) What regulator causes the transition between these attractors? How many variables change their activity between the initial and final state? what could be the biological meaning of this switch? how would you test this experimentally?

Exercise 8.13 (Biology) Finally, use the model to predict the rest of the attractor transitions because of single variable changes. Are all attractor transitions possible, or are there preferred differentiation routes? Why do you think this is?

8.5 References

• Rozum et al. (2023)
• Garcı́a-Gómez et al. (2020)

Garcı́a-Gómez, Mónica L, Diego Ornelas-Ayala, Adriana Garay-Arroyo, Berenice Garcı́a-Ponce, Marı́a de la Paz Sánchez, and Elena R Alvarez-Buylla. 2020. “A System-Level Mechanistic Explanation for Asymmetric Stem Cell Fates: Arabidopsis Thaliana Root Niche as a Study System.” Scientific Reports 10 (1): 3525.

Rozum, Jordan C, Colin Campbell, Eli Newby, Fatemeh Sadat Fatemi Nasrollahi, and Reka Albert. 2023. “Boolean Networks as Predictive Models of Emergent Biological Behaviors.” arXiv Preprint arXiv:2310.12901.