diff --git a/dingo/utils.py b/dingo/utils.py
index 473721bf..711531cf 100644
--- a/dingo/utils.py
+++ b/dingo/utils.py
@@ -204,13 +204,14 @@ def get_matrices_of_full_dim_polytope(A, b, Aeq, beq):
 
 
 
-def correlated_reactions(steady_states, pearson_cutoff = 0.90, indicator_cutoff = 10, 
+def correlated_reactions(steady_states, reactions=[], pearson_cutoff = 0.90, indicator_cutoff = 10, 
                          cells = 10, cop_coeff = 0.3, lower_triangle = True):
     """A Python function to calculate the pearson correlation matrix of a model
        and filter values based on the copula's indicator
 
     Keyword arguments:
-    steady_states -- A numpy array of the generated steady states fluxes 
+    steady_states -- A numpy array of the generated steady states fluxes
+    reactions -- A list with the model's reactions
     pearson_cutoff -- A cutoff to filter reactions based on pearson coefficient
     indicator_cutoff -- A cutoff to filter reactions based on indicator value
     cells -- Number of cells to compute the copula
@@ -256,6 +257,10 @@ def correlated_reactions(steady_states, pearson_cutoff = 0.90, indicator_cutoff
             np.fill_diagonal(filtered_corr_matrix, 1)
             return filtered_corr_matrix
     else:
+        # a dictionary that will store for each filtered reaction combination, 
+        # the pearson correlation value, the copula's indicator value 
+        # and the correlation classification
+        indicator_dict = {}
         
         # keep only the lower triangle
         corr_matrix = np.tril(corr_matrix)
@@ -264,17 +269,16 @@ def correlated_reactions(steady_states, pearson_cutoff = 0.90, indicator_cutoff
     
         # find indices of correlation matrix where correlation occurs
         corr_indices = np.argwhere((corr_matrix > pearson_cutoff) | (corr_matrix < -pearson_cutoff))
-            
-        # count reactions with positive or negative correlation based on indicator
-        positive = 0
-        negative = 0
     
         # compute copula for each set of correlated reactions
         for i in range(0, corr_indices.shape[0]):
             
             index1 = corr_indices[i][0]
             index2 = corr_indices[i][1]
-        
+            
+            reaction1 = reactions[index1]
+            reaction2 = reactions[index2]
+            
             flux1 = steady_states[index1]
             flux2 = steady_states[index2]
                 
@@ -302,26 +306,35 @@ def correlated_reactions(steady_states, pearson_cutoff = 0.90, indicator_cutoff
             # classify specific pair of reactions as positive or negative
             # correlated based on indicator cutoff 
             if indicator > indicator_cutoff:
-                positive += 1
+                pearson = filtered_corr_matrix[index1, index2]
+                indicator_dict[reaction1 + "_" + reaction2] = {'pearson': pearson,
+                                                               'indicator': indicator, 
+                                                               'classification': "positive"}
+
             elif indicator < 1/indicator_cutoff:
-                negative += 1
+                pearson = filtered_corr_matrix[index1, index2]
+                indicator_dict[reaction1 + "_" + reaction2] = {'pearson': pearson,
+                                                               'indicator': indicator, 
+                                                               'classification': "negative"}
+
             # if they do not overcome the cutoff replace their corresponding
             # value in the correlation matrix with 0
             else:
                 filtered_corr_matrix[index1, index2] = 0
                 filtered_corr_matrix[index2, index1] = 0
+                pearson = filtered_corr_matrix[index1, index2]
+                indicator_dict[reaction1 + "_" + reaction2] = {'pearson': pearson,
+                                                               'indicator': indicator, 
+                                                               'classification': "no correlation"}
                 
             print("Completed process of",i+1,"from",corr_indices.shape[0],"copulas")
-
-                
-        print(positive, "out of", i+1, "copulas were positive correlated based on copula indicator")
-        print(negative, "out of", i+1, "copulas were negative correlated based on copula indicator")
         
         if lower_triangle == True:
             filtered_corr_matrix[np.triu_indices(filtered_corr_matrix.shape[0], 1)] = np.nan
             np.fill_diagonal(filtered_corr_matrix, 1)
-            return filtered_corr_matrix
+            return filtered_corr_matrix, indicator_dict
 
         else:
             np.fill_diagonal(filtered_corr_matrix, 1)
-            return filtered_corr_matrix
\ No newline at end of file
+            return filtered_corr_matrix, indicator_dict
+            
\ No newline at end of file
diff --git a/tests/correlation.py b/tests/correlation.py
index 88239dd0..0e935bb0 100644
--- a/tests/correlation.py
+++ b/tests/correlation.py
@@ -14,11 +14,12 @@ def test_correlation(self):
         sampler = PolytopeSampler(dingo_model)
         steady_states = sampler.generate_steady_states()
         
-        # calculate correlation matrix without filtering from copula indicator
-        corr_matrix = correlated_reactions(steady_states, 
-                                           indicator_cutoff=2,
-                                           pearson_cutoff = 0.99999,
-                                           lower_triangle = False)
+        # calculate correlation matrix with filtering from copula indicator
+        corr_matrix, indicator_dict = correlated_reactions(steady_states,
+                                        reactions = reactions,
+                                        indicator_cutoff = 5,
+                                        pearson_cutoff = 0.999999,
+                                        lower_triangle = False)
 
         # sum values in the diagonal of the correlation matrix ==> 95*pearson ==> 95*1
         self.assertTrue(np.trace(corr_matrix) == len(reactions))
@@ -35,7 +36,7 @@ def test_correlation(self):
         
         # calculate correlation matrix without filtering from copula indicator
         corr_matrix = correlated_reactions(steady_states, 
-                                           indicator_cutoff=0,
+                                           indicator_cutoff = 0,
                                            pearson_cutoff = 0,
                                            lower_triangle = True)