workflow_LoopDetect_Matlab

% LoopDetect: Feedback Loop Detection in ODE models in MATLAB

0001 %% LoopDetect: Feedback Loop Detection in ODE models in MATLAB
0002 
0003 %%
0004 %
0005 % Copyright Katharina Baum, 2020.
0006 % Developed for MATLAB 2019b and higher.
0007 
0008 %% Installation
0009 %
0010 % Download and unzip the content of the folder 'LoopDetect_for_Matlab'.
0011 % Within the MATLAB session, navigate to this folder and work there or add the path of the
0012 % folder to MATLAB's search path for the current MATLAB session
0013 % by MATLAB's |addpath()| function. Deending on where you stored the folder,
0014 % its name could be something like '/Users/Desktop/LoopDetect' on Mac or
0015 % 'C:\matlab\LoopDetect' on Windows.
0016 
0017 % Retrieve the LoopDetect folder location when having navigated to the folder
0018 % within the MATLAB session.
0019 LoopDetect_Folder_Name = pwd;
0020 % Add this location to MATLAB's searchpath.
0021 addpath(LoopDetect_Folder_Name)
0022 
0023 %%
0024 % Attention: The LoopDetect content
0025 % will be added only for the current session. If restarting MATLAB, you
0026 % will have to perform this procedure again unless the LoopDetect folder
0027 % is stored in a location that already belongs to MATLAB's searchpath
0028 % (view it with MATLAB's |path()|
0029 % function), or you navigate to the folder and work within the folder itself.
0030 
0031 %% In brief and quick start
0032 %
0033 % The function suite LoopDetect enables determining all feedback loops of
0034 % an ordinary differential equation (ODE) system at user-defined values of
0035 % the model parameters and of the modelled variables.
0036 %
0037 % The following call reports (up to 10) feedback loops for an ODE system
0038 % determined by a function, here |func_POSm4()|, at variable values |s_star|
0039 % (here these are all equal to 1, column vector). The function values are
0040 % only allowed to depend on the variable values, other parameters and time
0041 % have to be fixed.
0042 %
0043 s_star=[1,1,1,1]';
0044 klin=ones(1,8); knonlin=[2.5,3];
0045 loop_list=find_loops_vset(@(x)func_POSm4(1,x,klin,knonlin),s_star,10);
0046 disp(loop_list{1})
0047 first_loop=loop_list{1}.loop{1}
0048 
0049 %%
0050 % Here, the first loop in |loop_list|, |first_loop|, is a negative feedback
0051 % loop (|sign| in the table: |-1|) of length 4
0052 % in that variable 4 regulates variable 1,
0053 % variable 1 regulates variable 2, variable 2 regulates variable 3, and
0054 % variable 3 regulates variable 4.
0055 
0056 %% Introduction
0057 %
0058 % Ordinary differential equation (ODE) models are used frequently to
0059 % mathematically represent biological systems. Feedback loops are important
0060 % regulatory features of biological systems and can give rise to different
0061 % dynamic behavior such as multistability for positive feedback loops or
0062 % oscillations for negative feedback loops.
0063 %
0064 % The feedback loops in an ODE system can be detected with the help of its
0065 % Jacobian matrix, the matrix of partial derivatives of the variables.
0066 % It captures all interactions between the variables and gives rise to the
0067 % interaction graph of the ODE model. In this graph, each modelled variable
0068 % is a node and non-zero entries in the Jacobian matrix are (weighted)
0069 % edges of the graph. Interactions can be positive or negative, according
0070 % to the sign of the Jacobian matrix entry.
0071 %
0072 % Directed path detection in this graph is used to determine all feedback
0073 % loops (in graphs also called cycles or circuits) of the system. They are
0074 % marked by a set of directed interactions forming a chain in which only
0075 % the first and the last node (variable) is the same. Thereby, self-loops
0076 % (loops of length one) can also occur.
0077 %
0078 % LoopDetect allows for detection of all loops of the graph and also reports
0079 % the sign of each loop, i.e. whether it is a positive feedback loop
0080 % (the number of negative interactions is even) or a negative feedback loop
0081 % (the number of negative interactions is uneven).
0082 % The output is a table that captures the order of the variables forming
0083 % the loop, their length and the sign of each loop.
0084 %
0085 % Except for solving the ODE system to generate variable values of interest
0086 % which requires the optimization toolbox, only MATLAB base functions are
0087 % employed and no further toolboxes are required. The algorithms rely on
0088 % Christopher Maes' implementation of
0089 % Johnson's algorithm for path detection (find_elem_circuits.m on github) and
0090 % on Brandon Kuczenski's implementation of Tarjan's algorithm for detecting
0091 % strongly connected components (tarjan.m, obtained from MATLAB File Exchange).
0092 %
0093 % Note that this tool was developed under MATLAB 2019b and that for earlier
0094 % MATLAB versions, it cannot be guaranteed that LoopDetect functions or parts
0095 % of this workflow run without errors.
0096 
0097 %% Solving the ODE model to generate variable values of interest
0098 
0099 klin=[165,0.044,0.27,550,5000,78,4.4,5.1];
0100 knonlin=[0.3,2];
0101 [t,sol]=ode15s(@(t,x)func_POSm4(t,x,klin,knonlin),[0,50],ones(1,4));
0102 s_star=sol(end,:); %the last point of the simulation is chosen
0103 
0104 %%
0105 % First, the ODE model is solved for a certain set of parameters, |klin|,
0106 % |knonlin|, with a
0107 % numerical solver in order to determine values of interest |s_star| for the
0108 % modelled variables of the ODE system. These could be steady state values,
0109 % values at a specific point in time (e.g. after a stimulus) or even a set
0110 % of parameter values (see section
0111 % Determining loops over a set of variable values). Thus, you can skip this step if
0112 % you already have a point of interest in state space, or if you want to
0113 % use dummy values such as |s_star=1:4|.
0114 
0115 %% Calculating the Jacobian matrix
0116 %
0117 % The supplied functions |numerical_jacobian()| and
0118 % |numerical_jacobian_complex()| can be used to determine numerically the Jacobian matrix
0119 % of an ODE system at a certain set of values for the variables, |s_star|.
0120 % The approach is that of finite differences (with real step) or
0121 % complex step approach, the latter of which is supposed to deliver more exact results [Martins et al., 2003] and
0122 % relies on the native implementation of complex numbers in MATLAB. The
0123 % input function handle |f| (in the example derived from |func_POSm4()|,
0124 % positive feedback chain model from [Baum et al., 2016]) points to a function that defines the time
0125 % derivatives of the modelled variables as a vector:
0126 % |f_i(s)=dS_i/dt|.
0127 % When supplied to |numerical_jacobian_complex()|, it is allowed to depend only on the values of
0128 % the variables; other dependencies, e.g. on parameter values, should be
0129 % removed by chosing fixed values.
0130 
0131 klin=[165,0.044,0.27,550,5000,78,4.4,5.1];
0132 knonlin=[0.3,2];
0133 j_matrix=numerical_jacobian_complex(@(s)func_POSm4(1,s,klin,knonlin),...
0134 s_star);
0135 signed_jacobian=sign(j_matrix)
0136     
0137 %%
0138 % The (i,j)th entry of the Jacobian matrix denotes the partial derivative
0139 % of variable |S_i| with respect to variable |S_j|,
0140 % |J_ij=delta S_i/delta S_j|,
0141 % which is positive if |S_j| has a direct positive effect on |S_i|, negative
0142 % if |S_j| has a direct negative effect on |S_i| and zero if |S_j| does not
0143 % have a direct effect on |S_i|. For example, the entry in row 1,
0144 % column 4, |J_14|, of the matrix above is |-1|, meaning that in the
0145 % underlying ODE model, variable 4 negatively regulates variable 1.
0146     
0147 %% Computing all feedback loops and useful functions for loop search
0148 %
0149 % The Jacobian matrix is used to compute feedback loops in the generated interaction graph.
0150 % For this, the default function is |find_loops()|, in that strongly
0151 % connected components are determined to reduce runtime. For smaller
0152 % systems, the function |find_loops_noscc()| skips this step and thus can be
0153 % faster. The optional second input argument sets an upper limit to the number of detected
0154 % and reported loops and thus can prevent overly long runtime (but also
0155 % potentially not all loops are returned).
0156 %
0157 
0158 loop_list=find_loops(j_matrix)
0159 
0160 %%
0161 % Single loops can be examined by entering the corresponding entry.
0162 
0163 for i=1:6
0164     disp(loop_list.loop{i})
0165 end
0166 %%
0167 % The function |loop_summary()| provides a convenient report on total number of
0168 % loops, subdivided by their lengths and signs.
0169 
0170 disp(loop_summary(loop_list))
0171 %%
0172 % One can filter the loop list for loops containing specific variables,
0173 % for example the one with index 2:
0174 noi = 2;
0175 loops_with_node2=loop_list(...
0176     cellfun(@(z) ismember(noi,z),loop_list.loop),:)
0177 
0178 %%
0179 % Search a loop list for loops containing specific edges defined by the
0180 % indices of the ingoing and outgoing nodes. This example returns the
0181 % indices of all loops with a regulation of node 3 by node 2. These are
0182 % only two here.
0183 %
0184 loop_edge_ind=find_edge(loop_list,2,3);
0185 loops_with_edge_2_to_3=loop_list(loop_edge_ind,:)
0186 for i=1:2
0187     disp(loops_with_edge_2_to_3.loop{i})
0188 end
0189 
0190 %%
0191 % Saving and reading loop lists from files can be done using MATLAB's
0192 % |save()| and |load()| functions. They keep the correct data format, but
0193 % objects have to be retrieved from a struct.
0194 
0195 save('loop_list_example.mat', 'loops_with_node2')
0196 loaded_loops_with_node2=load('loop_list_example.mat')
0197 % access the loaded data table from the struct (might not be required
0198 % in earlier MATLAB versions)
0199 loaded_loops_with_node2.loops_with_node2
0200 
0201 %%
0202 % Reading and writing loop lists to tabular format can be performed via
0203 % MATLAB's |writetable()| and |readtable()| functions. Here, we choose tabs as
0204 % delimiters. Note that the formatting is lost.
0205 
0206 writetable(loops_with_node2,'loop_list_example.txt','Delimiter','\t')
0207 loops_with_node2_readin = readtable('loop_list_example.txt')
0208 
0209 
0210 %% Computing feedback loops over multiple sets of variable values of interest
0211 %
0212 % In this example of a model of the bacterial cell cycle [Li et al., 2008],
0213 % we demonstrate how feedback loops can be determined over multiple sets of
0214 % variable values. Here, we focus on the solution of the ODE
0215 % systems along the time axis (provided in 'li08_solution.txt').
0216 
0217 % load sets of variable values (solution to the ODE over time)
0218 sol = readmatrix('li08_solution.txt');
0219 % possible alternative: sol = dlmread('li08_solution.txt');
0220 [loop_rep,loop_rep_index,jac_rep,jac_rep_index]=find_loops_vset(...
0221     @(x)func_li08(0,x), sol(:,2:end)', 1e5, false);
0222 
0223 %%
0224 % The solutions give rise to seven different loop lists.
0225 disp(length(loop_rep))
0226 
0227 %%
0228 % Here, two examples of resulting loop lists are given (without
0229 % self-loops).
0230 
0231 loop_list_2=loop_rep{2}(loop_rep{2}.length>1,:)
0232 loop_list_7=loop_rep{7}(loop_rep{7}.length>1,:)
0233 
0234 %%
0235 % These results could be plotted along the solution, e.g. by indicating a
0236 % certain background color for certain specific loop lists, and analyzed
0237 % further to discover reasons of changing loops. Please note that in order
0238 % to obtain the sample solution in |li08_solution.txt| also event functions
0239 % are required; the solution cannot be retrieved from integrating
0240 % |func_li08()| alone. Please refer to the model's publication [Li et al.,
0241 % 2008] for details.
0242 
0243 
0244 %% Comparing two loop lists
0245 %
0246 % We might want to compare loops of two systems, e.g. as obtained in the
0247 % example above. For this, we provide the function |compare_loop_list()|.
0248 % Thereby, loop indices in
0249 % the compared systems should point to the same variables, otherwise a meaningful
0250 % comparison is not possible. This could be the case for example if we
0251 % examine a system in which only one regulation has changed (for example
0252 % the positive feedback chain model vs. the negative feedback chain model, [Baum et al. 2016]),
0253 % or if regulations changed within one system when determining loops for
0254 % multiple sets of variables of interest (along a dynamic trajectory, at
0255 % different steady states of the system).
0256 
0257 % Function with positive regulation
0258 klin=[165,0.044,0.27,550,5000,78,4.4,5.1];
0259 knonlin=[0.3,2];
0260 j_matrix=numerical_jacobian_complex(@(s)func_POSm4(1,s,klin,knonlin),...
0261     s_star);
0262 loop_list_pos=find_loops(j_matrix);
0263 
0264 % Function with negative regulation. The altered regulation affects
0265 % two entries of the Jacobian matrix. Parameter values and the set of
0266 % variable values remain identical.
0267 j_matrix_neg=j_matrix;
0268 j_matrix_neg(1:2,4)=-j_matrix(1:2,4);
0269 loop_list_neg=find_loops(j_matrix_neg);
0270 
0271 % compute comparison
0272 [ind_a_id,ind_a_switch,ind_a_notin,ind_b_id,ind_b_switch]=...
0273     compare_loop_list(loop_list_pos,loop_list_neg);
0274 
0275 %%
0276 % Only the four self-loops remain identical in both systems.
0277 
0278 disp(loop_list_pos(ind_a_id,:))
0279 for i=1:4
0280     disp(loop_list_pos.loop{ind_a_id(i)})
0281 end
0282 
0283 %%
0284 % These are the corresponding loops in the negatively regulated system.
0285 
0286 disp(loop_list_neg(ind_b_id,:))
0287 for i=1:4
0288     disp(loop_list_neg.loop{ind_b_id(i)})
0289 end
0290 
0291 %%
0292 % Two loops are the same in both systems but they have switched their signs.
0293 
0294 loops_switch_pos=loop_list_pos(ind_a_switch,:)
0295 for i=1:2
0296     disp(loops_switch_pos.loop{i})
0297 end
0298 loops_switch_neg=loop_list_neg(ind_b_switch,:)
0299 for i=1:2
0300     disp(loops_switch_neg.loop{i})
0301 end
0302 
0303 %%
0304 % All loops in the positively regulated system do also occur in the negatively regulated
0305 % system, i.e. |ind_a_notin| is empty.
0306 
0307 ind_a_notin
0308 
0309 
0310 %% List of functions in LoopDetect - Alphabetical overview
0311 %
0312 % * |compare_loop_list| - compare two loop lists
0313 % * |find_edge| - find loops with a certain edge in a loop list
0314 % * |find_loops_noscc| - detect feedback loops from a Jacobian or adjacency
0315 % matrix without determining strongly connected components, suitable for
0316 % smaller or densely connected systems
0317 % * |find_loops_vset| - detect feedback loops from a function and sets of
0318 % values for the modelled variables
0319 % * |find_loops| - detect feedback loops from a Jacobian or adjacency matrix
0320 % * |loop_summary| - summarize loops list contents by lengths and signs
0321 % * |numerical_jacobian_complex| - finite difference numerical determination of the Jacobian
0322 % matrix with a complex step, provides higher accuracy than
0323 % |numerical_jacobian|
0324 % * |numerical_jacobian| - finite difference numerical determination of the Jacobian
0325 % matrix
0326 % * |sort_loop_index| - sort the reported loops to always start with the
0327 % lowest node index
0328 %
0329 
0330 %% References
0331 %
0332 % Baum K, Politi AZ, Kofahl B, Steuer R, Wolf J. Feedback, Mass
0333 % Conservation and Reaction Kinetics Impact the Robustness of Cellular
0334 % Oscillations. PLoS Comput Biol. 2016;12(12):e1005298.
0335 %
0336 % Brandon Kuczenski (2018). tarjan(e)
0337 % (https://www.mathworks.com/matlabcentral/fileexchange/50707-tarjan-e),
0338 % MATLAB Central File Exchange. Retrieved August 14, 2018.
0339 %
0340 % Li S, Brazhnik P, Sobral B, Tyson JJ. A Quantitative Study of the
0341 % Division Cycle of Caulobacter crescentus Stalked Cells. Plos Comput Biol.
0342 % 2008;4(1):e9.
0343 %
0344 % Christopher Maes (2011). find_elem_circuits.m
0345 % https://gist.github.com/cmaes/1260153
0346 %
0347 % Martins JRRA, Sturdza P, Alonso JJ. The complex-step derivative
0348 % approximation. ACM Trans Math Softw. 2003;29(3):245–62.
0349 %