diff --git a/src/libaevol/7/ExpManager_7.cpp b/src/libaevol/7/ExpManager_7.cpp
index 9427f061b140912bedcee86710eeaab55195be20..15630a4162d44bb617bc9f5f8aaf166b59fe600b 100644
--- a/src/libaevol/7/ExpManager_7.cpp
+++ b/src/libaevol/7/ExpManager_7.cpp
@@ -90,6 +90,11 @@ ExpManager_7::ExpManager_7(ExpManager* exp_m) {
     dna_size[indiv_id] = current_individuals[indiv_id]->dna_->length();
   }
 
+#ifdef __REGUL
+  phenotypic_target_handler_ = new SIMD_PhenotypicTargetHandler_R(
+        dynamic_cast<PhenotypicTargetHandler_R*>(exp_m->world()->phenotypic_target_handler()),
+        exp_m->exp_s());
+#else
 #ifdef PHENOTYPE_VECTOR
   target = new double[PHENOTYPE_VECTOR_SIZE];
     for (int i = 0; i < PHENOTYPE_VECTOR_SIZE; i++) {
@@ -99,6 +104,7 @@ ExpManager_7::ExpManager_7(ExpManager* exp_m) {
 #else
   target = new Vector_Fuzzy(*(Fuzzy*)(exp_m->world()->phenotypic_target_handler()->phenotypic_target().fuzzy()));
 #endif
+#endif
 
 #ifdef WITH_PERF_TRACES
   std::ofstream perf_traces_file_;
@@ -132,28 +138,112 @@ void ExpManager_7::selection(int indiv_id) {
 
   int cur_x,cur_y;
 
+#ifdef __REGUL
+  double ** fitness_sum_local_tab_;
+#endif
+
+  if (fitness_function == FITNESS_LOCAL_SUM) {
+#ifdef __REGUL
+    fitness_sum_local_tab_ = new double*[fitness_function_scope_x*fitness_function_scope_y];
+        for (int tab_id = 0; tab_id < fitness_function_scope_x*fitness_function_scope_y; tab_id++)
+          fitness_sum_local_tab_[tab_id] = new double[phenotypic_target_handler_->nb_env_];
+
+        for (int env_id = 0; env_id < phenotypic_target_handler_->nb_env_; env_id++) {
+
+          int tab_id = 0;
+          fitness_sum_local_tab_[tab_id][env_id] = 0;
+
+          for (int8_t i = -1; i < selection_scope_x - 1; i++) {
+            for (int8_t j = -1; j < selection_scope_y - 1; j++) {
+              cur_x = (x + i + grid_width) % grid_width;
+              cur_y = (y + j + grid_height) % grid_height;
+
+              int16_t new_x,new_y;
+              for (int8_t ii = -1; ii < fitness_function_scope_x - 1; ii++) {
+                for (int8_t jj = -1; jj < fitness_function_scope_y - 1; jj++) {
+                  //TODO: Check values HERE !
+
+                  new_x = (cur_x + ii + grid_width) % grid_width;
+                  new_y = (cur_y + jj + grid_height) % grid_height;
+
+                  fitness_sum_local_tab_[tab_id][env_id] +=
+                      prev_internal_simd_struct[new_x * grid_height + new_y]->fitness_by_env_id_[env_id];
+                }
+              }
+              tab_id++;
+            }
+          }
+        }
+#else
+    printf("Fitness local sum is not supported for Aevol (only R-Aevol)\n");
+    exit(-1);
+#endif
+  }
+
+  int tab_id = 0;
+
   for (int8_t i = -1 ; i < selection_scope_x-1 ; i++) {
     for (int8_t j = -1; j < selection_scope_y - 1; j++) {
       cur_x = (x + i + grid_width) % grid_width;
       cur_y = (y + j + grid_height) % grid_height;
 
-
-      local_fit_array[count] =
-          previous_individuals[cur_x * grid_height + cur_y]->fitness;
-
-      local_meta_array[count] =
-          previous_individuals[cur_x * grid_height + cur_y]->metaerror;
-
-      if (local_meta_array[count] !=
-          previous_individuals[cur_x * grid_height + cur_y]->metaerror) {
-        printf("NONONNONONONN\n");
+      if (fitness_function == FITNESS_EXP) {
+        local_fit_array[count] =
+            prev_internal_simd_struct[cur_x * grid_height + cur_y]->fitness;
+        local_meta_array[count] =
+            prev_internal_simd_struct[cur_x * grid_height + cur_y]
+                ->metaerror;
+      } else if (fitness_function == FITNESS_GLOBAL_SUM) {
+#ifdef __REGUL
+        double composed_fitness = 0;
+            for (int env_id = 0; env_id < phenotypic_target_handler_->nb_env_;
+                 env_id++) {
+              composed_fitness +=
+                  prev_internal_simd_struct[cur_x * grid_height + cur_y]->fitness_by_env_id_[env_id] /
+                  fitness_sum_tab_[env_id];
+            }
+            composed_fitness /= phenotypic_target_handler_->nb_env_;
+            local_fit_array[count] = composed_fitness;
+#else
+        printf("Fitness local sum is not supported for Aevol (only R-Aevol)\n");
+        exit(-1);
+#endif
+      } else if (fitness_function == FITNESS_LOCAL_SUM) {
+#ifdef __REGUL
+        double composed_fitness = 0;
+            for (int env_id = 0; env_id < phenotypic_target_handler_->nb_env_;
+                 env_id++) {
+              composed_fitness +=
+                  prev_internal_simd_struct[cur_x * grid_height + cur_y]->fitness_by_env_id_[env_id] /
+                  fitness_sum_local_tab_[tab_id][env_id];
+            }
+            composed_fitness /= phenotypic_target_handler_->nb_env_;
+            local_fit_array[count] = composed_fitness;
+#else
+        printf("Fitness local sum is not supported for Aevol (only R-Aevol)\n");
         exit(-1);
+#endif
       }
+
       sum_local_fit += local_fit_array[count];
+
       count++;
+      tab_id++;
     }
   }
 
+
+  if (fitness_function == FITNESS_LOCAL_SUM) {
+#ifdef __REGUL
+    for (int tab_id = 0; tab_id < fitness_function_scope_x * fitness_function_scope_y; tab_id++)
+          delete[] fitness_sum_local_tab_[tab_id];
+        delete[] fitness_sum_local_tab_;
+#else
+    printf("Fitness local sum is not supported for Aevol (only R-Aevol)\n");
+    exit(-1);
+#endif
+  }
+
   for(int16_t i = 0 ; i < neighborhood_size ; i++) {
     probs[i] = local_fit_array[i]/sum_local_fit;
   }
@@ -444,7 +534,10 @@ ExpManager_7::~ExpManager_7() {
 
   delete dna_factory_;
 
+#ifdef __REGUL
+#else
   delete target;
+#endif
 }
 
 
@@ -1216,6 +1309,10 @@ void ExpManager_7::compute_protein(int indiv_id) {
       }
     }
   }
+
+#ifdef __REGUL
+  ((SIMD_List_Metadata*)internal_simd_struct[indiv_id]->metadata_)->add_inherited_proteins();
+#endif
 }
 
 
@@ -1748,6 +1845,22 @@ void ExpManager_7::build_phenotypic_target(PhenotypicTargetHandler* phenotypic_t
 }
 
 void ExpManager_7::compute_fitness(int indiv_id, double selection_pressure) {
+#ifdef __REGUL
+  Vector_Fuzzy* delta = new Vector_Fuzzy(*internal_simd_struct[indiv_id]->phenotype);
+      delta->sub(phenotypic_target_handler_->targets_fuzzy_[env_id]);
+
+      if (phenotypic_target_handler_->var_method_ == SWITCH_IN_A_LIST)
+        internal_simd_struct[indiv_id]->metaerror_by_env_id_[env_id] = delta->get_geometric_area();
+      else if (phenotypic_target_handler_->var_method_ == ONE_AFTER_ANOTHER)
+        internal_simd_struct[indiv_id]->metaerror_by_env_id_[env_id] += delta->get_geometric_area();
+
+      delete delta;
+
+      if (phenotypic_target_handler_->var_method_ == SWITCH_IN_A_LIST)
+        internal_simd_struct[indiv_id]->fitness_by_env_id_[env_id] = exp(
+            -selection_pressure *
+            ((double) internal_simd_struct[indiv_id]->metaerror_by_env_id_[env_id]));
+#else
 #ifdef PHENOTYPE_VECTOR
   for (int fuzzy_idx = 0; fuzzy_idx < PHENOTYPE_VECTOR_SIZE; fuzzy_idx++) {
 
@@ -1780,15 +1893,267 @@ void ExpManager_7::compute_fitness(int indiv_id, double selection_pressure) {
         current_individuals[indiv_id]->fitness = exp(
       -selection_pressure *
       ((double)current_individuals[indiv_id]->metaerror));
+#endif
 }
 
 
+#ifdef __REGUL
+void ExpManager_7::compute_network(int indiv_id, double selection_pressure) {
+  // Allocate fitness and metarror array
+  if (phenotypic_target_handler_->var_method_ == SWITCH_IN_A_LIST) {
+    internal_simd_struct[indiv_id]->fitness_by_env_id_ = new double[phenotypic_target_handler_->nb_eval_];
+    internal_simd_struct[indiv_id]->metaerror_by_env_id_ = new double[phenotypic_target_handler_->nb_eval_];
+  } else if (phenotypic_target_handler_->var_method_ == ONE_AFTER_ANOTHER) {
+    internal_simd_struct[indiv_id]->fitness_by_env_id_ = new double[phenotypic_target_handler_->nb_env_];
+    internal_simd_struct[indiv_id]->metaerror_by_env_id_ = new double[phenotypic_target_handler_->nb_env_];
+
+    for (int env_id = 0; env_id < phenotypic_target_handler_->nb_env_; env_id++) {
+      internal_simd_struct[indiv_id]->fitness_by_env_id_[env_id] = 0;
+      internal_simd_struct[indiv_id]->metaerror_by_env_id_[env_id] = 0;
+    }
+  }
+
+  // Add signals
+  ((SIMD_List_Metadata*)internal_simd_struct[indiv_id]->metadata_)->signal_proteins_.resize(phenotypic_target_handler_->signals_models_.size());
+  int i = 0;
+  for (auto signal : phenotypic_target_handler_->signals_models_) {
+    int glob_prot_idx = internal_simd_struct[indiv_id]->metadata_->proteins_count();
+    internal_simd_struct[indiv_id]->metadata_->set_proteins_count(
+        internal_simd_struct[indiv_id]->metadata_->proteins_count() +
+        1);
+
+    pProtein* prot = new pProtein(signal);
+    internal_simd_struct[indiv_id]->metadata_->protein_add(glob_prot_idx, prot);
+    ((SIMD_List_Metadata*)internal_simd_struct[indiv_id]->metadata_)->signal_proteins_[i] = (prot);
+    i++;
+  }
+
+  // Set influences
+  internal_simd_struct[indiv_id]->metadata_->rna_begin();
+  for (int i = 0; i < internal_simd_struct[indiv_id]->metadata_->rna_count(); i++) {
+    pRNA* rna = internal_simd_struct[indiv_id]->metadata_->rna_next();
+    if (rna != nullptr) {
+      if (rna->is_coding_) {
+        rna->nb_influences_ = 0;
+
+        int32_t enhancer_position = rna->enhancer_position(internal_simd_struct[indiv_id]->dna_->length());
+        int32_t operator_position = rna->operator_position(internal_simd_struct[indiv_id]->dna_->length());
+
+        double enhance,operate;
+        internal_simd_struct[indiv_id]->metadata_->protein_begin();
+        for (int j = 0; j < internal_simd_struct[indiv_id]->metadata_->proteins_count(); j++) {
+          pProtein* prot =
+              internal_simd_struct[indiv_id]->metadata_->protein_next();
+          if (prot != nullptr) {
+            enhance = rna->affinity_with_protein(enhancer_position, prot, internal_simd_struct[indiv_id], exp_m_);
+            operate = rna->affinity_with_protein(operator_position, prot, internal_simd_struct[indiv_id], exp_m_);
+
+            if (enhance != 0.0 || operate != 0.0) {
+
+              rna->affinity_list.push_back(AffinityFactor(prot->e,enhance,operate));
+              prot->is_TF_ = true;
+
+              rna->nb_influences_ ++;
+            }
+          }
+        }
+      }
+    }
+  }
+}
+
+void ExpManager_7::update_network(int indiv_id, double selection_pressure) {
+
+  internal_simd_struct[indiv_id]->metadata_->rna_begin();
+  for (int x = 0; x < internal_simd_struct[indiv_id]->metadata_->rna_count(); x++) {
+    pRNA* rna = internal_simd_struct[indiv_id]->metadata_->rna_next();
+    if (rna != nullptr) {
+      if (rna->is_coding_) {
+
+        double enhancer_activity = 0;
+        double operator_activity = 0;
+
+        for (auto affinity: rna->affinity_list) {
+          enhancer_activity +=
+              affinity.enhancer_factor * affinity.protein_concentration;
+          operator_activity +=
+              affinity.operator_factor * affinity.protein_concentration;
+        }
+
+        ProteinConcentration enhancer_activity_pow_n =
+            enhancer_activity == 0
+                ? 0
+                : pow(enhancer_activity, exp_m_->exp_s()->get_hill_shape_n());
+        ProteinConcentration operator_activity_pow_n =
+            operator_activity == 0
+                ? 0
+                : pow(operator_activity, exp_m_->exp_s()->get_hill_shape_n());
+        rna->synthesis_rate =
+            rna->e *
+            (exp_m_->exp_s()->get_hill_shape() /
+             (operator_activity_pow_n + exp_m_->exp_s()->get_hill_shape())) *
+            (1 + ((1 / rna->e) - 1) * (enhancer_activity_pow_n /
+                                       (enhancer_activity_pow_n +
+                                        exp_m_->exp_s()->get_hill_shape())));
+      }
+    }
+  }
+
+  internal_simd_struct[indiv_id]->metadata_->protein_begin();
+  for (int j = 0; j < internal_simd_struct[indiv_id]->metadata_->proteins_count(); j++) {
+    pProtein* prot =
+        internal_simd_struct[indiv_id]->metadata_->protein_next();
+    if (!prot->signal_) {
+      prot->delta_concentration_ = 0;
+
+      for (auto rna : prot->rna_list_) {
+        prot->delta_concentration_ += rna->synthesis_rate;
+      }
+
+      prot->delta_concentration_ -= exp_m_->exp_s()->get_degradation_rate() * prot->e;
+      prot->delta_concentration_ *= 1.0/exp_m_->exp_s()->get_nb_degradation_step();
+    }
+  }
+
+  // Apply the changes in concentrations we have just computed
+
+  internal_simd_struct[indiv_id]->metadata_->protein_begin();
+  for (int j = 0; j < internal_simd_struct[indiv_id]->metadata_->proteins_count(); j++) {
+    pProtein* prot =
+        internal_simd_struct[indiv_id]->metadata_->protein_next();
+    if (!prot->signal_) prot->e += prot->delta_concentration_;
+  }
+}
+
+void ExpManager_7::evaluate_network(int indiv_id, double selection_pressure, int env_id) {
+  update_phenotype(indiv_id);
+
+  compute_fitness(indiv_id,selection_pressure,env_id);
+
+  if (phenotypic_target_handler_->var_method_ == SWITCH_IN_A_LIST)
+    internal_simd_struct[indiv_id]->metaerror += internal_simd_struct[indiv_id]->metaerror_by_env_id_[env_id];
+}
+
+void ExpManager_7::finalize_network(int indiv_id, double selection_pressure) {
+  if (phenotypic_target_handler_->var_method_ == SWITCH_IN_A_LIST) {
+    internal_simd_struct[indiv_id]->metaerror =
+        internal_simd_struct[indiv_id]->metaerror /
+        (double)(phenotypic_target_handler_->nb_indiv_age_);
+  }
+  else if (phenotypic_target_handler_->var_method_ == ONE_AFTER_ANOTHER) {
+    for (int env_id = 0; env_id < phenotypic_target_handler_->nb_env_; env_id++) {
+      internal_simd_struct[indiv_id]->metaerror_by_env_id_[env_id] =
+          internal_simd_struct[indiv_id]->metaerror_by_env_id_[env_id] / 10.0;
+      internal_simd_struct[indiv_id]->metaerror +=
+          internal_simd_struct[indiv_id]->metaerror_by_env_id_[env_id];
+
+      internal_simd_struct[indiv_id]->fitness_by_env_id_[env_id] = exp(
+          -selection_pressure *
+          ((double) internal_simd_struct[indiv_id]->metaerror_by_env_id_[env_id]));
+    }
+    internal_simd_struct[indiv_id]->metaerror =
+        internal_simd_struct[indiv_id]->metaerror / (double) (phenotypic_target_handler_->nb_env_);
+  }
+
+  internal_simd_struct[indiv_id]->fitness = exp(
+      -selection_pressure *
+      ((double) internal_simd_struct[indiv_id]->metaerror));
+}
+
+void ExpManager_7::solve_network(int indiv_id, double selection_pressure) {
+  internal_simd_struct[indiv_id]->metadata_->protein_begin();
+  for (int j = 0;
+       j < internal_simd_struct[indiv_id]->metadata_->proteins_count(); j++) {
+    pProtein* prot =
+        internal_simd_struct[indiv_id]->metadata_->protein_next();
+    if (!prot->signal_) {
+      prot->e = prot->initial_e_;
+    }
+  }
+
+  compute_network(indiv_id, selection_pressure);
+
+  internal_simd_struct[indiv_id]->metaerror = 0;
+
+  if (phenotypic_target_handler_->var_method_ == ONE_AFTER_ANOTHER) {
+    for (int16_t env_i = 0; env_i < phenotypic_target_handler_->nb_env_; env_i++) {
+
+      //Set the concentration of signals for this age
+      for (auto prot1: ((SIMD_List_Metadata*)internal_simd_struct[indiv_id]->metadata_)->signal_proteins_) {
+        prot1->e = 0;
+      }
+
+      for (auto prot_id : phenotypic_target_handler_->env_signals_list_[env_i]) {
+        ((SIMD_List_Metadata*)internal_simd_struct[indiv_id]->metadata_)->signal_proteins_[prot_id]->e = 0.9;
+      }
+
+      for (int16_t i = 0; i < 10; i++) {
+
+        for (int j = 0; j < 10; j++) {
+          update_network(indiv_id,selection_pressure);
+        }
+
+        // If we have to evaluate the individual at this age
+        evaluate_network(indiv_id,selection_pressure,env_i);
+      }
+    }
+
+    finalize_network(indiv_id,selection_pressure);
+  } else {
+    std::set<int>* eval = exp_m_->exp_s()->get_list_eval_step();
+    // i is thus the age of the individual
+    for (int16_t i = 1; i <= phenotypic_target_handler_->nb_indiv_age_; i++) {
+      //Set the concentration of signals for this age
+      for (auto prot1: ((SIMD_List_Metadata*)internal_simd_struct[indiv_id]->metadata_)->signal_proteins_) {
+        prot1->e = 0;
+      }
+
+      for (auto prot_id : phenotypic_target_handler_->env_signals_list_[phenotypic_target_handler_->list_env_id_[i]]) {
+        ((SIMD_List_Metadata*)internal_simd_struct[indiv_id]->metadata_)->signal_proteins_[prot_id]->e = 0.9;
+      }
+
+      for (int j = 0; j < exp_m_->exp_s()->get_nb_degradation_step(); j++) {
+        update_network(indiv_id,selection_pressure);
+      }
+
+      // If we have to evaluate the individual at this age
+      if (eval->find(i) != eval->end()) {
+        evaluate_network(indiv_id,selection_pressure, phenotypic_target_handler_->list_env_id_[i]);
+      }
+    }
+
+    finalize_network(indiv_id,selection_pressure);
+  }
+}
+
+void ExpManager_7::update_phenotype( int indiv_id ) {
+  delete internal_simd_struct[indiv_id]->phenotype;
+  compute_phenotype(indiv_id);
+}
+#endif
+
+
 void ExpManager_7::run_a_step(double w_max, double selection_pressure,bool optim_prom) {
 #pragma omp single
   {
     nb_clones_ = 0;
     cumulate_size = 0;
     cumulate_diff = 0;
+
+
+    if (exp_m_->sel()->fitness_func() == FITNESS_GLOBAL_SUM) {
+#ifdef __REGUL
+      fitness_sum_tab_ = new double[phenotypic_target_handler_->nb_env_];
+        for (int env_id = 0; env_id < phenotypic_target_handler_->nb_env_; env_id++) {
+          fitness_sum_tab_[env_id] = 0;
+          for (int indiv_id = 0; indiv_id < exp_m_->nb_indivs(); indiv_id++)
+            fitness_sum_tab_[env_id] += prev_internal_simd_struct[indiv_id]->fitness_by_env_id_[env_id];
+        }
+#else
+      printf("Fitness local sum is not supported for Aevol (only R-Aevol)\n");
+      exit(-1);
+#endif
+    }
   }
 
 #pragma omp for schedule(dynamic)
diff --git a/src/libaevol/7/ExpManager_7.h b/src/libaevol/7/ExpManager_7.h
index e79945c08f836548daa3b371a67459b709aa388c..5d9b5f16b8c11fddd1dffd2969215ef4ff9e70d5 100644
--- a/src/libaevol/7/ExpManager_7.h
+++ b/src/libaevol/7/ExpManager_7.h
@@ -80,7 +80,14 @@ class ExpManager_7 : public Observable{
 
 
   void build_phenotypic_target(PhenotypicTargetHandler* phenotypic_target_handler);
-
+#ifdef __REGUL
+  void compute_network(int indiv_id, double selection_pressure);
+  void update_network(int indiv_id, double selection_pressure);
+  void evaluate_network(int indiv_id, double selection_pressure, int env_id);
+  void finalize_network(int indiv_id, double selection_pressure);
+  void solve_network(int indiv_id, double selection_pressure);
+  void update_phenotype( int indiv_id );
+#endif
 
   void set_stats(Stats* stats) { stats_ = stats; }
 
@@ -103,10 +110,15 @@ class ExpManager_7 : public Observable{
  private:
   ExpManager* exp_m_;
   int* dna_size;
+#ifdef __REGUL
+  Vector_Fuzzy** targets;
+SIMD_PhenotypicTargetHandler_R* phenotypic_target_handler_;
+#else
 #ifdef PHENOTYPE_VECTOR
   double* target;
 #else
   Vector_Fuzzy* target;
+#endif
 #endif
   bool first_gener_ = true;
 
diff --git a/src/libaevol/7/List_Metadata.cpp b/src/libaevol/7/List_Metadata.cpp
index fb6190b4a3f2df569c4383156efd7fa89539c4ad..f8c3667736fdd6f7fa3e89e6c0d6d381a7def611 100644
--- a/src/libaevol/7/List_Metadata.cpp
+++ b/src/libaevol/7/List_Metadata.cpp
@@ -934,4 +934,4 @@ namespace aevol {
     void List_Metadata::proteins_clear() {
         proteins_.clear();
     }
-}
\ No newline at end of file
+}
diff --git a/src/libaevol/7/List_Metadata.h b/src/libaevol/7/List_Metadata.h
index 147145dbf30d7b5542e46a7bae0f8562ece6fd38..08b16bc77e94ed20e71b18f26371c70832c3b761 100644
--- a/src/libaevol/7/List_Metadata.h
+++ b/src/libaevol/7/List_Metadata.h
@@ -73,6 +73,8 @@ namespace aevol {
             }
         };
 
+        void add_inherited_proteins();
+
         /** Getter **/
 
         void set_iterators() {
diff --git a/src/libaevol/7/Protein_7.cpp b/src/libaevol/7/Protein_7.cpp
index ab9fb754d82903e24101c69c7b74cb4d80fad2eb..227e78626e0d054fd95488c236b35f6d419c4fb2 100644
--- a/src/libaevol/7/Protein_7.cpp
+++ b/src/libaevol/7/Protein_7.cpp
@@ -26,3 +26,32 @@
 
 
 #include "Protein_7.h"
+Protein_7::Protein_7(Protein_R* prot_sig) {
+  protein_length = prot_sig->length();
+  e = prot_sig->concentration();
+
+  signal_ = prot_sig->is_signal();
+
+  nb_codons_ = prot_sig->AA_list().size();
+
+  for (int i = 0; i < nb_codons_; i++)
+    codon_list[i] = prot_sig->_cod_tab[i];
+}
+
+
+Protein_7::Protein_7(pProtein* prot) {
+  protein_start   = prot->protein_start;
+  protein_end     = prot->protein_end;
+  protein_length  = prot->protein_length;
+  leading_lagging = prot->leading_lagging;
+  e               = prot->e;
+  is_init_        = true;
+  signal_         = prot->signal_;
+  nb_codons_      = prot->nb_codons_;
+
+  for (int i = 0; i < nb_codons_; i++)
+    codon_list[i] = prot->codon_list[i];
+
+  initial_e_ = prot->initial_e_;
+  inherited_ = true;
+}
diff --git a/src/libaevol/7/Protein_7.h b/src/libaevol/7/Protein_7.h
index 968de1278b95a4f74eb282375a5068d546c4952c..5ff84e3a8af639c1b14a2c9041493ba5e38015c3 100644
--- a/src/libaevol/7/Protein_7.h
+++ b/src/libaevol/7/Protein_7.h
@@ -31,6 +31,7 @@
 #include <cstdint>
 
 namespace aevol {
+
 class Protein_7 {
  public:
   Protein_7(){};
@@ -38,15 +39,27 @@ class Protein_7 {
             int32_t t_protein_end,
             int32_t t_protein_length,
             int8_t t_leading_lagging,
-            double t_e) {
+            double t_e, Rna_7* rna) {
     protein_start   = t_protein_start;
     protein_end     = t_protein_end;
     protein_length  = t_protein_length;
     leading_lagging = t_leading_lagging;
     e               = t_e;
     is_init_        = true;
+
+#ifdef __REGUL
+    rna_list_.push_back(rna);
+      initial_e_ = e;
+#endif
   }
 
+#ifdef __REGUL
+  Protein_7(Protein_R* prot_sig);
+
+
+  Protein_7(pProtein* prot);
+#endif
+
   bool operator<(const Protein_7& other);
 
   int32_t protein_start;
@@ -60,6 +73,15 @@ class Protein_7 {
   bool is_functional;
 
   bool is_init_ = false;
+#ifdef __REGUL
+  bool is_TF_;
+
+  double initial_e_ = -1;
+  double    delta_concentration_;
+  bool      inherited_ = false;
+  bool      signal_;
+  std::list<pRNA*> rna_list_;
+#endif
 };
 }
 
diff --git a/src/libaevol/7/Rna_7.cpp b/src/libaevol/7/Rna_7.cpp
index 215f84f492f436ee7df4bd67edf8f7a964018ea2..77bc1302b38b35f1c2c145fc8b40c5fa1f846748 100644
--- a/src/libaevol/7/Rna_7.cpp
+++ b/src/libaevol/7/Rna_7.cpp
@@ -26,3 +26,51 @@
 
 
 #include "Rna_7.h"
+
+
+double Rna_7::affinity_with_protein( int32_t index, pProtein *protein, Internal_SIMD_Struct* indiv, ExpManager* exp_m ) {
+  int32_t len = protein->protein_length;
+
+  if (len > 5) {
+    double max = 0;
+    double temp = 1;
+
+    int32_t quadon_tab[5];
+
+    for (int32_t pos = index; pos < index+5; pos++) {
+
+      int8_t quadon[4];
+
+      if (leading_lagging == LEADING) {
+        for (int8_t i = 0; i < QUADON_SIZE; i++) {
+          quadon[i] = (indiv->dna_->get_lead(pos + i) == '1')
+                      ? 1 << (QUADON_SIZE - i - 1)
+                      : 0;
+        }
+      } else {
+        for (int8_t i = 0; i < QUADON_SIZE; i++) {
+          quadon[i] = (indiv->dna_->get_lag(pos - i) != '1')
+                      ? (QUADON_SIZE - i - 1)
+                      : 0;
+        }
+      }
+
+      quadon_tab[pos - index] = quadon[0] + quadon[1] + quadon[2] + quadon[3];
+    }
+
+    for (int32_t i = 0; i < len - 4; i++) {
+      temp = 1;
+
+      for (int8_t j = 0; j < 5; j++) {
+        temp *= exp_m->exp_s()->get_binding_matrix(quadon_tab[j],
+                                                   protein->codon_list[i + j]);
+      }
+
+      max = (max < temp) ? temp : max;
+    }
+
+    return max;
+  } else {
+    return 0.0;
+  }
+}
diff --git a/src/libaevol/7/Rna_7.h b/src/libaevol/7/Rna_7.h
index 56451898e9eb338b998561b06b545ea7893e7604..a9045982d2db30f82532aa99764c93398d900851 100644
--- a/src/libaevol/7/Rna_7.h
+++ b/src/libaevol/7/Rna_7.h
@@ -33,6 +33,21 @@
 #include <cstdint>
 
 namespace aevol {
+
+class AffinityFactor {
+ public:
+  AffinityFactor(double pconcentration, double efactor, double ofactor) {
+    protein_concentration = pconcentration;
+    enhancer_factor = efactor;
+    operator_factor = ofactor;
+  }
+
+  double protein_concentration;
+  double enhancer_factor;
+  double operator_factor;
+};
+
+
 class Rna_7 {
  public:
   Rna_7(){};
@@ -53,6 +68,50 @@ class Rna_7 {
 
   ~Rna_7() {}
 
+
+#ifdef __REGUL
+  std::list<AffinityFactor> affinity_list;
+
+    int nb_influences_ = 0;
+
+    int32_t enhancer_position(int32_t length) {
+      if(leading_lagging == LEADING)
+      {
+        return (begin - 20)  % ( length ) < 0 ?
+               ((begin - 20)  % ( length )) + ( length ) :
+               (begin - 20)  % ( length );
+      }
+      else  // strand_ = LAGGING
+      {
+        return (begin + 20)  % ( length ) < 0 ?
+               ((begin + 20)  % ( length )) + ( length ) :
+               (begin + 20)  % ( length );
+      }
+    }
+
+  int32_t operator_position(int32_t length) {
+    if(leading_lagging == LEADING)
+    {
+      return (begin + PROM_SIZE)  % ( length ) < 0 ?
+             (begin + PROM_SIZE)  % ( length ) + (length) :
+             (begin + PROM_SIZE)  % ( length );
+    }
+    else  // strand_ = LAGGING
+    {
+      return (begin - PROM_SIZE)  % ( length ) < 0 ?
+             (begin - PROM_SIZE)  % ( length ) + (length) :
+             (begin - PROM_SIZE)  % ( length );
+    }
+  }
+
+  double affinity_with_protein( int32_t index, pProtein *protein,
+                                     Internal_SIMD_Struct* indiv,
+                               ExpManager* exp_m);
+
+  double synthesis_rate;
+#endif
+
+
   int32_t begin;
   int32_t end;
   int8_t leading_lagging; // 0 = leading, 1 = lagging
diff --git a/src/libaevol/CMakeLists.txt b/src/libaevol/CMakeLists.txt
index bc82eaad1867ea0c1871ee399af2d9553ab1ab42..b770e8691457fcd0a4f1067c0675d9644eb5b612 100644
--- a/src/libaevol/CMakeLists.txt
+++ b/src/libaevol/CMakeLists.txt
@@ -202,9 +202,9 @@ if ( ${with-raevol} )
           raevol/Protein_R.h
           raevol/Protein_R.cpp
           raevol/Rna_R.h
-          raevol/Rna_R.cpp
-	  #"SIMD_Abstract Metadata.cpp" "SIMD_Abstract Metadata.h" Map_Metadata.cpp Map_Metadata.h
-            7/ExpManager_7.cpp 7/ExpManager_7.h 7/Promoter.cpp 7/Promoter.h 7/Rna_7.cpp 7/Rna_7.h 7/Protein_7.cpp 7/Protein_7.h)
+          raevol/Rna_R.cpp 
+	  #"SIMD_Abstract Metadata.cpp" "SIMD_Abstract Metadata.h" SIMD_Map_Metadata.cpp SIMD_Map_Metadata.h
+            raevol/SIMD_PhenotypicTargetHandler_R.cpp raevol/SIMD_PhenotypicTargetHandler_R.h)
 endif ()
 
 add_library(aevol ${aevol_sources})
diff --git a/src/libaevol/population/Selection.cpp b/src/libaevol/population/Selection.cpp
index c40ff6c0efef2edd8454856b679f4d26afed1744..8e8c1be94115ed21608b5de5d1cc39c25060ae92 100644
--- a/src/libaevol/population/Selection.cpp
+++ b/src/libaevol/population/Selection.cpp
@@ -944,14 +944,6 @@ if (exp_m_->record_tree() || exp_m_->light_tree()) {
   return new_indiv;
 }
 
-void Selection::run_life(Individual* new_indiv) {
-  // Evaluate new individual
-  new_indiv->Evaluate();
-
-  // Compute statistics
-  new_indiv->compute_statistical_data();
-
-}
 
 #ifdef __REGUL
 void Selection::run_life(Individual_R* new_indiv) {
@@ -967,6 +959,15 @@ void Selection::run_life(Individual_R* new_indiv) {
     // Compute statistics
     new_indiv->compute_statistical_data();
 
+}
+#else
+void Selection::run_life(Individual* new_indiv) {
+  // Evaluate new individual
+  new_indiv->Evaluate();
+
+  // Compute statistics
+  new_indiv->compute_statistical_data();
+
 }
 #endif
 
diff --git a/src/libaevol/raevol/PhenotypicTargetHandler_R.cpp b/src/libaevol/raevol/PhenotypicTargetHandler_R.cpp
index 657eac8038387c4a8797e0f10fe5bf1cc7388bec..be577b8fb30a44ea0deeba91c7830948ea369d32 100644
--- a/src/libaevol/raevol/PhenotypicTargetHandler_R.cpp
+++ b/src/libaevol/raevol/PhenotypicTargetHandler_R.cpp
@@ -474,36 +474,6 @@ void PhenotypicTargetHandler_R::BuildPhenotypicTargetModel( int16_t id) {
       }
       phenotypic_target->fuzzy()->add_point(new_point.x, new_point.y);
     }
-
-    if (FuzzyFactory::fuzzyFactory->get_fuzzy_flavor() == 1) {
-      HybridFuzzy* fuz = (HybridFuzzy*) phenotypic_target->fuzzy();
-
-      for (int i = 1; i < fuz->get_pheno_size(); i++) {
-        if (fuz->points()[i] == 0.0) {
-          int minL = i - 1;
-          int maxL = i + 1;
-          int dist = 1;
-
-          while (fuz->points()[maxL] == 0.0) {
-            maxL++;
-            dist++;
-          }
-          double inc = 0.0;
-          if (fuz->points()[maxL] > fuz->points()[minL]) {
-            inc = (fuz->points()[maxL] - fuz->points()[minL]) / dist;
-          } else {
-            inc = (fuz->points()[minL] - fuz->points()[maxL]) / dist;
-            minL = maxL;
-          }
-
-          for (int j = i; j < maxL; j++) {
-            fuz->points()[j] = fuz->points()[minL] + inc;
-            inc += inc;
-          }
-
-        }
-      }
-    }
   }
   // Add lower and upper bounds
   phenotypic_target->fuzzy()->clip(AbstractFuzzy::min, Y_MIN);
diff --git a/src/libaevol/raevol/PhenotypicTargetHandler_R.h b/src/libaevol/raevol/PhenotypicTargetHandler_R.h
index 8459b47fec652459702d6f13cbb74a922ad97a6f..9b2238ba992abf87c28524a33f6f7ce223068122 100644
--- a/src/libaevol/raevol/PhenotypicTargetHandler_R.h
+++ b/src/libaevol/raevol/PhenotypicTargetHandler_R.h
@@ -179,6 +179,14 @@ class PhenotypicTargetHandler_R : public virtual PhenotypicTargetHandler
 
     void ShuffleRandomlySignals();
 
+  std::vector<PhenotypicTarget_R*> phenotypic_targets_;
+  std::vector<std::list<Gaussian>> env_gaussians_list_;
+  std::vector<std::list<int16_t>> env_signals_list_;
+  std::vector<Protein_R*> signals_models_;
+  std::list<Protein_R*> signals_models_list_;
+  double env_switch_probability_;
+  int16_t _nb_indiv_age;
+
  protected :
   // ==========================================================================
   //                            Protected Methods
@@ -197,13 +205,6 @@ class PhenotypicTargetHandler_R : public virtual PhenotypicTargetHandler
   //                               Attributes
   // ==========================================================================
 
-  std::vector<PhenotypicTarget_R*> phenotypic_targets_;
-  std::vector<std::list<Gaussian>> env_gaussians_list_;
-  std::vector<std::list<int16_t>> env_signals_list_;
-  std::vector<Protein_R*> signals_models_;
-  std::list<Protein_R*> signals_models_list_;
-  double env_switch_probability_;
-  int16_t _nb_indiv_age;
 
   bool hasChanged_;
 
diff --git a/src/libaevol/raevol/SIMD_PhenotypicTargetHandler_R.cpp b/src/libaevol/raevol/SIMD_PhenotypicTargetHandler_R.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2dfd21d55d2b2f1a75af5c0e3d311d0256b89acc
--- /dev/null
+++ b/src/libaevol/raevol/SIMD_PhenotypicTargetHandler_R.cpp
@@ -0,0 +1,132 @@
+//
+// Created by arrouan on 18/08/2020.
+//
+
+#include "SIMD_PhenotypicTargetHandler_R.h"
+#include "ExpSetup.h"
+namespace aevol {
+
+SIMD_PhenotypicTargetHandler_R::SIMD_PhenotypicTargetHandler_R(PhenotypicTargetHandler_R *handler, ExpSetup* exp_s) {
+  env_gaussians_list_.resize(handler->env_gaussians_list_.size());
+  nb_env_ = env_gaussians_list_.size();
+
+  int i = 0;
+  for (auto gaussian_env : handler->env_gaussians_list_) {
+    for (auto gauss : gaussian_env) {
+      env_gaussians_list_[i].push_back(Gaussian(gauss));
+      i++;
+    }
+  }
+
+  env_signals_list_.resize(handler->env_signals_list_.size());
+  i=0;
+  for (auto env_signal : handler->env_signals_list_) {
+    for (auto e_signal : env_signal) {
+      env_signals_list_[i].push_back(e_signal);
+      i++;
+    }
+  }
+
+  i = 0;
+  for (auto signal_protein : handler->signals_models_) {
+    pProtein* prot = new pProtein(signal_protein);
+    signals_models_[i] = prot;
+    i++;
+  }
+
+  env_switch_probability_ = handler->env_switch_probability_;
+  nb_indiv_age_ = handler->_nb_indiv_age;
+
+  var_method_ = handler->var_method();
+  var_prng_ = handler->var_prng_;
+
+  sampling_ = handler->sampling();
+
+  targets_fuzzy_by_id_ = new Vector_Fuzzy*[nb_env_];
+
+  for (int env_id = 0; env_id < nb_env_; env_id++) {
+    targets_fuzzy_by_id_[env_id] = new Vector_Fuzzy();
+
+    if (not env_gaussians_list_.at(env_id).empty()) {
+      for (int16_t i = 0; i <= sampling_; i++) {
+        Point new_point = Point(
+            X_MIN + (double) i * (X_MAX - X_MIN) / (double) sampling_, 0.0);
+        int gi = 0;
+        for (const Gaussian& g: env_gaussians_list_.at(env_id)) {
+          gi++;
+          new_point.y += g.compute_y(new_point.x);
+        }
+        targets_fuzzy_by_id_[env_id]->add_point(new_point.x, new_point.y);
+      }
+    }
+    // Add lower and upper bounds
+    targets_fuzzy_by_id_[env_id]->clip(AbstractFuzzy::min, Y_MIN);
+    targets_fuzzy_by_id_[env_id]->clip(AbstractFuzzy::max, Y_MAX);
+
+    // Simplify (get rid of useless points)
+    targets_fuzzy_by_id_[env_id]->simplify();
+  }
+
+  if (var_method_ == SWITCH_IN_A_LIST) {
+    targets_fuzzy_ = new Vector_Fuzzy*[nb_indiv_age_];
+
+    if (nb_env_ <= 1) {
+      for (int age = 0; age < nb_indiv_age_; age++) {
+        targets_fuzzy_[age] = targets_fuzzy_by_id_[0];
+      }
+    }
+
+    nb_eval_ = exp_s->get_list_eval_step()->size();
+  }
+
+
+}
+
+void SIMD_PhenotypicTargetHandler_R::ApplyVariation() {
+
+  switch (var_method_) {
+  case NO_VAR :
+    return;
+  case SWITCH_IN_A_LIST : {
+    if ( nb_env_ <= 1 ) {
+      break;
+    }
+
+    int16_t* list_of_old_target_id = list_env_id_;
+    list_env_id_ = new int16_t[nb_indiv_age_];
+
+    // Shortcuts used
+    int16_t id_new_env = list_of_old_target_id[nb_indiv_age_-1];
+    int16_t id_old_env = list_of_old_target_id[nb_indiv_age_-1];
+
+    hasChanged_ = false;
+
+    for (int16_t i = 0; i < nb_indiv_age_ ; i++) {
+        // if we have to change of environment :
+        double env_chang = var_prng_->random();
+
+        if (env_chang < env_switch_probability_) {
+          //we have to change to a new env that have an id different from the old one
+          while (id_new_env == id_old_env) {
+            id_new_env = var_prng_->random(nb_env_);
+          }
+          //The environment has changed
+          id_old_env = id_new_env;
+        }
+
+        list_env_id_[i] = id_new_env;
+        targets_fuzzy_[i] = targets_fuzzy_by_id_[id_new_env];
+
+        if (list_env_id_[i] != list_of_old_target_id[i])
+          hasChanged_ = true;
+    }
+    break; }
+  case ONE_AFTER_ANOTHER:
+    break;
+  default :
+    Utils::ExitWithDevMsg("Unknown variation method", __FILE__, __LINE__);
+    break;
+  }
+}
+
+}
\ No newline at end of file
diff --git a/src/libaevol/raevol/SIMD_PhenotypicTargetHandler_R.h b/src/libaevol/raevol/SIMD_PhenotypicTargetHandler_R.h
new file mode 100644
index 0000000000000000000000000000000000000000..469968c67ce39e22fea099248f773dd3de077456
--- /dev/null
+++ b/src/libaevol/raevol/SIMD_PhenotypicTargetHandler_R.h
@@ -0,0 +1,54 @@
+//
+// Created by arrouan on 18/08/2020.
+//
+
+#ifndef AEVOL_SIMD_PHENOTYPICTARGETHANDLER_R_H
+#define AEVOL_SIMD_PHENOTYPICTARGETHANDLER_R_H
+
+#include "population/SIMD_Individual.h"
+#include "PhenotypicTargetHandler_R.h"
+#include "phenotype/Gaussian.h"
+#include "raevol/Protein_R.h"
+#include "Vector_Fuzzy.h"
+
+#include <list>
+#include <vector>
+
+namespace aevol {
+class pProtein;
+
+class SIMD_PhenotypicTargetHandler_R {
+ public:
+  SIMD_PhenotypicTargetHandler_R(PhenotypicTargetHandler_R* handler, ExpSetup* exp_s);
+
+  void ApplyVariation();
+
+  std::vector<pProtein*> signals_models_;
+  std::vector<std::list<int16_t>> env_signals_list_;
+
+  Vector_Fuzzy** targets_fuzzy_;
+  PhenotypicTargetVariationMethod var_method_;
+
+  int16_t* list_env_id_;
+
+  int16_t nb_indiv_age_;
+  int16_t nb_eval_;
+  int16_t nb_env_;
+ protected:
+  std::vector<std::list<Gaussian>> env_gaussians_list_;
+
+
+  double env_switch_probability_;
+  int16_t sampling_;
+
+
+  Vector_Fuzzy** targets_fuzzy_by_id_;
+
+  std::shared_ptr<JumpingMT> var_prng_;
+
+  bool hasChanged_ = false;
+};
+
+}
+
+#endif //AEVOL_SIMD_PHENOTYPICTARGETHANDLER_R_H