Reseach Projects- Inverse Perturbation for Optimal Intervention in Gene Regulatory Networks



Research Projects




Inverse Perturbation for Optimal Intervention in Gene Regulatory Networks

Analysis and intervention in the dynamics of gene regulatory networks is at the heart of emerging efforts in the development of modern treatment of numerous ailments including cancer. The ultimate goal is to develop methods to intervene in the function of living organisms in order to drive cells away from a malignant state into a benign form. A serious limitation of much of the previous work in cancer network analysis is the use of external control, which requires intervention at each time step, for an indefinite time interval. This is in sharp contrast to the proposed approach, which relies on the solution of an inverse perturbation problem to introduce a one-time intervention in the structure of regulatory networks. This isolated intervention transforms the attractor set of the dynamic system to a unique attractor characterized by the desired steady-state distribution.

We formulate the optimal intervention problem in gene regulatory networks as a minimal-perturbation of the network in order to force it to converge to a desired-steady-state distribution of gene regulation. We cast optimal intervention in gene regulation as a convex optimization problem, thus providing a globally optimal solution which can be efficiently computed for networks with thousands of nodes. The criteria adopted for optimality is chosen to minimize potential adverse effects as a consequence of the intervention strategy. We consider a perturbation that minimizes (i) the overall energy of change between the original and controlled networks and (ii) the time needed to reach the desired steady-state of gene regulation. Furthermore, we show that there is an inherent tradeoff between minimizing the energy of the perturbation and the convergence rate to the desired distribution. We apply the proposed control to the Human melanoma gene regulatory network.


DATA: Download the probability transition matrix of the Human melanoma network:

CODE: The MATLAB code uses CVX, a package for specifying and solving convex programs. One needs to download and install the CVX package before running the MATLAB program
Min_energy_perturbation.m . For CVX installation instructions, please follow the guidelines in

Download the MATLAB code:
Min_energy_perturbation.m with output the minimal-perturbation energy matrix.xls

Download the log file of the code:
log file.txt

The code is also displayed below:
0001   %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0002   % The program "Min_energy_perturbation.m" reads a square matrix, the 
0003   % probability transition matrix of the network, "P_matrix.xls" in excel format and returns the outputs "C", "nu", 
0004   % "gamma", and "mu" defined as follows: 
0005   % C: the minimal energy perturbation matrix, which forces the network to 
0006   % transition from its initial steady-state to the desired steady-state pid. 
0007   % nu: the Lagrange multiplier associated with the equality: pid'*(P+C) == 
0008   % pid' 
0009   % gamma: the Lagrange multiplier associated with the equality C*OneVec == 
0010   % ZeroVec 
0011   % mu: the Lagrange multiplier associated with the inequality P+C >= 0. 
0012   % 
0013   % This code uses the CVX software. For installation instructions, please 
0014   % refer to: 
0015   %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0018   %% Read the probability transition matrix "P_matrix.xls" 
0019   [P,Txt,Raw]=xlsread('P_matrix');
0021   %% Global variable definition 
0022   N = size(P,1);
0023   OneVec = ones(N,1);
0024   ZeroVec = zeros(N,1);
0026   %% Define the desired steady-state \pi_d 
0027   pid1 = 0.015525 * ones(N/2,1);
0028   pid2 = 10^(-4) * ones(N/2,1);
0029   pid = [pid1;pid2];
0031   %% The optimization algorithm using the CVX software 
0032   cvx_begin
0033       variable C(N,N);
0034       dual variable nu;
0035       dual variable gamma;
0036       dual variable mu;
0038       minimize( norm(C) );
0039       subject to
0040       nu: pid'*(P+C) == pid'; 
0041       gamma: C*OneVec == ZeroVec; 
0042       mu: P+C >= 0; 
0043   cvx_end

Download the MATLAB code for Fig. 1: plot_SLEM.m

Download the MATLAB code for Fig. 2:

Download the MATHEMATICA code for Fig. 3: