Chemicals interaction and degradation

Hi everyone, I’m trying to implement the interaction between different chemicals in a BDM model. More specifically, the model should be able to replicate the depletion of a chemical (the extracellular matrix, ECM) mediated by another one (the matrix metalloproteinases, MMP), with the amount of depletion dependent on the concentration of the two chemicals and on a binding coefficient between the two.
What I’ve tried so far is the following: at every time step a Standalone Operation goes through every box in the diffusion grid of the MMP, gets the concentration of ECM and MMP in that box, then computes the amount of ECM depleted by the MMP and finally reduces the concentration of the ECM in that box by this amount.
This is the code that I’m using:

struct Depletion : public StandaloneOperationImpl {
  BDM_OP_HEADER(Depletion);

  void operator()() override {
    auto* sim = Simulation::GetActive();
    auto* rm = sim->GetResourceManager();

    auto* depletesDg = rm->GetDiffusionGrid(depletes_);
    auto* toBeDepletedDg = rm->GetDiffusionGrid(toBeDepleted_);
    if (depletesDg->GetResolution() >= toBeDepletedDg->GetResolution()) {
      auto boxes_depletes = depletesDg->GetNumBoxesArray();
      auto box_length_depletes = depletesDg->GetBoxLength();
      auto grid_dim_depletes = depletesDg->GetDimensions();
      auto box_volume_depletes = depletesDg->GetBoxVolume();
      auto box_volume_toBeDepleted = toBeDepletedDg->GetBoxVolume();
      double volume_ratio = box_volume_depletes / box_volume_toBeDepleted;

      auto nx = boxes_depletes[0];
      auto ny = boxes_depletes[1];
      auto nz = boxes_depletes[2];

#pragma omp parallel for
      for (size_t x = 0; x < nx; x++) {
        double real_x = grid_dim_depletes[0] + x * box_length_depletes;
        for (size_t y = 0; y < ny; y++) {
          double real_y = grid_dim_depletes[2] + y * box_length_depletes;
          for (size_t z = 0; z < nz; z++) {
            double real_z = grid_dim_depletes[4] + z * box_length_depletes;
            Double3 depletionPosition = {real_x, real_y, real_z};
            auto depletesConc = depletesDg->GetConcentration(depletionPosition);
            auto toBeDepletedConc =
                toBeDepletedDg->GetConcentration(depletionPosition);

            if (depletesConc > 0 && toBeDepletedConc > 0) {
              double amount = depletesConc * toBeDepletedConc * binding_rate_ *
                              volume_ratio;

              if (amount > toBeDepletedConc) 
              {
                amount = toBeDepletedConc;
              }
              toBeDepletedDg->ChangeConcentrationBy(depletionPosition, -amount);
            }
          }
        }
      }
    } else {
      Log::Fatal(
          "Resolution of the depleted diffusion grid is bigger than that of "
          "the depletes diffusion grid");
    }
  }

  Substances toBeDepleted_, depletes_;
  double binding_rate_;
};

At the moment the code above isn’t giving the expected results: am I missing something? Is there a better way of implementing such a mechanism?

Thank you,

Nicolò

UPDATE: adding the depletion term (MMP*binding coefficient, where MMP is the average concentration of MMP) to the decay constant and then using this sum to set the decay constant of the chemical which is depleted (i.e. the ECM, via ModelInitializer::DefineSubstance(…)) at the beginning of the simulation leads to results which are way closer to the expected ones.