diff --git a/include/backend/dd/DDSimDebug.hpp b/include/backend/dd/DDSimDebug.hpp
index 3696a5e..777b1f6 100644
--- a/include/backend/dd/DDSimDebug.hpp
+++ b/include/backend/dd/DDSimDebug.hpp
@@ -57,9 +57,10 @@ struct DDSimulationState {
   std::vector<size_t> callReturnStack;
   std::map<size_t, std::map<std::string, std::string>> callSubstitutions;
   std::vector<std::pair<size_t, size_t>> restoreCallReturnStack;
-  std::map<size_t, std::vector<size_t>> dataDependencies;
+  std::map<size_t, std::vector<std::pair<size_t, size_t>>> dataDependencies;
+  std::map<size_t, std::set<size_t>> functionCallers;
   std::set<size_t> breakpoints;
-  std::vector<std::set<std::string>> targetQubits;
+  std::vector<std::vector<std::string>> targetQubits;
 
   bool paused;
 
@@ -127,8 +128,14 @@ bool checkAssertion(DDSimulationState* ddsim,
                     std::unique_ptr<Assertion>& assertion);
 std::string getClassicalBitName(DDSimulationState* ddsim, size_t index);
 size_t variableToQubit(DDSimulationState* ddsim, const std::string& variable);
+std::pair<size_t, size_t> variableToQubitAt(DDSimulationState* ddsim,
+                                            const std::string& variable,
+                                            size_t instruction);
 bool isSubStateVectorLegal(const Statevector& full,
                            std::vector<size_t>& targetQubits);
 std::vector<std::vector<Complex>>
 getPartialTraceFromStateVector(const Statevector& sv,
                                const std::vector<size_t>& traceOut);
+
+std::vector<std::string> getTargetVariables(DDSimulationState* ddsim,
+                                            size_t instruction);
diff --git a/include/backend/dd/DDSimDiagnostics.hpp b/include/backend/dd/DDSimDiagnostics.hpp
index f7ae200..4be6b59 100644
--- a/include/backend/dd/DDSimDiagnostics.hpp
+++ b/include/backend/dd/DDSimDiagnostics.hpp
@@ -8,6 +8,7 @@
 #include <map>
 #include <memory>
 #include <set>
+#include <vector>
 
 struct DDSimulationState;
 
@@ -17,12 +18,15 @@ struct DDDiagnostics {
   DDSimulationState* simulationState;
   std::map<size_t, std::set<size_t>> zeroControls;
   std::map<size_t, std::set<size_t>> nonZeroControls;
+
+  std::map<size_t, std::set<std::vector<size_t>>> actualQubits;
 };
 
 size_t dddiagnosticsGetNumQubits(Diagnostics* self);
 size_t dddiagnosticsGetInstructionCount(Diagnostics* self);
-Result dddiagnosticsInit([[maybe_unused]] Diagnostics* self);
+Result dddiagnosticsInit(Diagnostics* self);
 Result dddiagnosticsGetDataDependencies(Diagnostics* self, size_t instruction,
+                                        bool includeCallers,
                                         bool* instructions);
 Result dddiagnosticsGetInteractions(Diagnostics* self, size_t beforeInstruction,
                                     size_t qubit, bool* qubitsAreInteracting);
diff --git a/include/backend/diagnostics.h b/include/backend/diagnostics.h
index 5bef392..0685234 100644
--- a/include/backend/diagnostics.h
+++ b/include/backend/diagnostics.h
@@ -24,7 +24,7 @@ struct Diagnostics {
   size_t (*getNumQubits)(Diagnostics* self);
   size_t (*getInstructionCount)(Diagnostics* self);
   Result (*getDataDependencies)(Diagnostics* self, size_t instruction,
-                                bool* instructions);
+                                bool includeCallers, bool* instructions);
   Result (*getInteractions)(Diagnostics* self, size_t beforeInstruction,
                             size_t qubit, bool* qubitsAreInteracting);
   Result (*getZeroControlInstructions)(Diagnostics* self, bool* instructions);
diff --git a/include/common/parsing/CodePreprocessing.hpp b/include/common/parsing/CodePreprocessing.hpp
index 0855d71..b3d84db 100644
--- a/include/common/parsing/CodePreprocessing.hpp
+++ b/include/common/parsing/CodePreprocessing.hpp
@@ -5,8 +5,8 @@
 #include <cstddef>
 #include <map>
 #include <memory>
-#include <set>
 #include <string>
+#include <utility>
 #include <vector>
 
 struct Block {
@@ -18,7 +18,7 @@ struct Instruction {
   size_t lineNumber;
   std::string code;
   std::unique_ptr<Assertion> assertion;
-  std::set<std::string> targets;
+  std::vector<std::string> targets;
 
   size_t originalCodeStartPosition;
   size_t originalCodeEndPosition;
@@ -32,13 +32,13 @@ struct Instruction {
 
   std::map<std::string, std::string> callSubstitution;
 
-  std::vector<size_t> dataDependencies;
+  std::vector<std::pair<size_t, size_t>> dataDependencies;
 
   Block block;
   std::vector<size_t> childInstructions;
   Instruction(size_t inputLineNumber, std::string inputCode,
               std::unique_ptr<Assertion>& inputAssertion,
-              std::set<std::string> inputTargets, size_t startPos,
+              std::vector<std::string> inputTargets, size_t startPos,
               size_t endPos, size_t successor, bool isFuncCall,
               std::string function, bool inFuncDef, bool isFuncDef,
               Block inputBlock);
diff --git a/include/common/parsing/Utils.hpp b/include/common/parsing/Utils.hpp
index 830ec8f..f0552e4 100644
--- a/include/common/parsing/Utils.hpp
+++ b/include/common/parsing/Utils.hpp
@@ -14,3 +14,5 @@ std::string replaceString(std::string str, const std::string& from,
                           const std::string& to);
 
 std::string removeWhitespace(std::string str);
+
+bool variablesEqual(const std::string& v1, const std::string& v2);
diff --git a/src/backend/dd/DDSimDebug.cpp b/src/backend/dd/DDSimDebug.cpp
index 52e5b90..f90bd2b 100644
--- a/src/backend/dd/DDSimDebug.cpp
+++ b/src/backend/dd/DDSimDebug.cpp
@@ -24,6 +24,7 @@
 #include <cstring>
 #include <exception>
 #include <iostream>
+#include <iterator>
 #include <memory>
 #include <numeric>
 #include <random>
@@ -854,6 +855,52 @@ Result destroyDDSimulationState(DDSimulationState* self) {
 }
 
 //-----------------------------------------------------------------------------------------
+
+std::vector<std::string> getTargetVariables(DDSimulationState* ddsim,
+                                            size_t instruction) {
+  std::vector<std::string> result;
+  size_t parentFunction = -1ULL;
+  size_t i = instruction;
+  while (true) {
+    if (ddsim->functionDefinitions.find(i) !=
+        ddsim->functionDefinitions.end()) {
+      parentFunction = i;
+      break;
+    }
+    if (ddsim->instructionTypes[i] == RETURN) {
+      break;
+    }
+    if (i == 0) {
+      break;
+    }
+    i--;
+  }
+
+  const auto parameters = parentFunction != -1ULL
+                              ? ddsim->targetQubits[parentFunction]
+                              : std::vector<std::string>{};
+  for (const auto& target : ddsim->targetQubits[instruction]) {
+    if (std::find(parameters.begin(), parameters.end(), target) !=
+        parameters.end()) {
+      result.push_back(target);
+      continue;
+    }
+    const auto foundRegister =
+        std::find_if(ddsim->qubitRegisters.begin(), ddsim->qubitRegisters.end(),
+                     [target](const QubitRegisterDefinition& reg) {
+                       return reg.name == target;
+                     });
+    if (foundRegister != ddsim->qubitRegisters.end()) {
+      for (size_t j = 0; j < foundRegister->size; j++) {
+        result.push_back(target + "[" + std::to_string(j) + "]");
+      }
+    } else {
+      result.push_back(target);
+    }
+  }
+  return result;
+}
+
 size_t variableToQubit(DDSimulationState* ddsim, const std::string& variable) {
   auto declaration = replaceString(variable, " ", "");
   declaration = replaceString(declaration, "\t", "");
@@ -892,6 +939,42 @@ size_t variableToQubit(DDSimulationState* ddsim, const std::string& variable) {
   throw std::runtime_error("Unknown variable name " + var);
 }
 
+std::pair<size_t, size_t> variableToQubitAt(DDSimulationState* ddsim,
+                                            const std::string& variable,
+                                            size_t instruction) {
+  size_t sweep = instruction;
+  size_t functionDef = -1ULL;
+  while (sweep < ddsim->instructionTypes.size()) {
+    if (std::find(ddsim->functionDefinitions.begin(),
+                  ddsim->functionDefinitions.end(),
+                  sweep) != ddsim->functionDefinitions.end()) {
+      functionDef = sweep;
+      break;
+    }
+    if (ddsim->instructionTypes[sweep] == RETURN) {
+      break;
+    }
+    sweep--;
+  }
+
+  if (functionDef == -1ULL) {
+    // In the global scope, we can just use the register's index.
+    return {variableToQubit(ddsim, variable), functionDef};
+  }
+
+  // In a gate-local scope, we have to define qubit indices relative to the
+  // gate.
+  const auto& targets = ddsim->targetQubits[functionDef];
+
+  const auto found = std::find(targets.begin(), targets.end(), variable);
+  if (found == targets.end()) {
+    throw std::runtime_error("Unknown variable name " + variable);
+  }
+
+  return {static_cast<size_t>(std::distance(targets.begin(), found)),
+          functionDef};
+}
+
 double complexMagnitude(Complex& c) {
   return std::sqrt(c.real * c.real + c.imaginary * c.imaginary);
 }
@@ -1283,6 +1366,7 @@ std::string preprocessAssertionCode(const char* code,
   ddsim->qubitRegisters.clear();
   ddsim->successorInstructions.clear();
   ddsim->dataDependencies.clear();
+  ddsim->functionCallers.clear();
   ddsim->targetQubits.clear();
 
   for (auto& instruction : instructions) {
@@ -1293,7 +1377,17 @@ std::string preprocessAssertionCode(const char* code,
     ddsim->instructionEnds.push_back(instruction.originalCodeEndPosition);
     ddsim->dataDependencies.insert({instruction.lineNumber, {}});
     for (const auto& dependency : instruction.dataDependencies) {
-      ddsim->dataDependencies[instruction.lineNumber].push_back(dependency);
+      ddsim->dataDependencies[instruction.lineNumber].emplace_back(
+          dependency.first, dependency.second);
+    }
+    if (instruction.isFunctionCall) {
+      const size_t successorInFunction = instruction.successorIndex;
+      const size_t functionIndex = successorInFunction - 1;
+      if (ddsim->functionCallers.find(functionIndex) ==
+          ddsim->functionCallers.end()) {
+        ddsim->functionCallers.insert({functionIndex, {}});
+      }
+      ddsim->functionCallers[functionIndex].insert(instruction.lineNumber);
     }
 
     // what exactly we do with each instruction depends on its type:
diff --git a/src/backend/dd/DDSimDiagnostics.cpp b/src/backend/dd/DDSimDiagnostics.cpp
index af1edbc..a89d1d8 100644
--- a/src/backend/dd/DDSimDiagnostics.cpp
+++ b/src/backend/dd/DDSimDiagnostics.cpp
@@ -9,6 +9,7 @@
 #include <algorithm>
 #include <cstddef>
 #include <cstdint>
+#include <iterator>
 #include <map>
 #include <memory>
 #include <set>
@@ -59,9 +60,90 @@ size_t dddiagnosticsGetInstructionCount(Diagnostics* self) {
       &ddd->simulationState->interface);
 }
 
-Result dddiagnosticsInit([[maybe_unused]] Diagnostics* self) { return OK; }
+Result dddiagnosticsInit(Diagnostics* self) {
+  auto* ddd = toDDDiagnostics(self);
+  ddd->zeroControls.clear();
+  ddd->nonZeroControls.clear();
+  ddd->actualQubits.clear();
+  return OK;
+}
+
+size_t findReturn(DDSimulationState* state, size_t instruction) {
+  size_t current = instruction;
+  while (state->instructionTypes[current] != RETURN) {
+    current++;
+  }
+  return current;
+}
+
+void visitCall(DDSimulationState* ddsim, size_t current, size_t qubitIndex,
+               std::set<size_t>& visited, std::set<size_t>& toVisit) {
+  const auto gateStart = ddsim->successorInstructions[current];
+  const auto gateDefinition = gateStart - 1;
+  const std::string stringToSearch =
+      ddsim->targetQubits[gateDefinition][qubitIndex];
+  auto checkInstruction = findReturn(ddsim, gateStart);
+  while (checkInstruction >= gateStart) {
+    const auto found =
+        std::find(ddsim->targetQubits[checkInstruction].begin(),
+                  ddsim->targetQubits[checkInstruction].end(), stringToSearch);
+    if (ddsim->instructionTypes[checkInstruction] != RETURN &&
+        found != ddsim->targetQubits[checkInstruction].end()) {
+      if (visited.find(checkInstruction) == visited.end()) {
+        toVisit.insert(checkInstruction);
+      }
+      if (ddsim->instructionTypes[checkInstruction] == CALL) {
+        const auto position =
+            std::distance(ddsim->targetQubits[checkInstruction].begin(), found);
+        visitCall(ddsim, checkInstruction, static_cast<size_t>(position),
+                  visited, toVisit);
+      }
+      break;
+    }
+    if (checkInstruction == 0) {
+      break;
+    }
+    checkInstruction--;
+  }
+}
+
+std::set<size_t> getUnknownCallers(DDSimulationState* ddsim,
+                                   size_t instruction) {
+  std::set<size_t> unknownCallers;
+
+  std::set<size_t> toVisit{};
+  std::set<size_t> visited{};
+
+  while (true) {
+    instruction--;
+    if (ddsim->functionDefinitions.find(instruction) !=
+        ddsim->functionDefinitions.end()) {
+      unknownCallers.insert(instruction);
+      for (const auto caller : ddsim->functionCallers[instruction]) {
+        if (visited.find(caller) == visited.end()) {
+          toVisit.insert(caller);
+        }
+      }
+    }
+
+    if (instruction == 0 || ddsim->instructionTypes[instruction] == RETURN ||
+        ddsim->functionDefinitions.find(instruction) !=
+            ddsim->functionDefinitions.end()) {
+      if (toVisit.empty()) {
+        break;
+      }
+
+      instruction = *toVisit.begin();
+      toVisit.erase(instruction);
+      visited.insert(instruction);
+    }
+  }
+
+  return unknownCallers;
+}
 
 Result dddiagnosticsGetDataDependencies(Diagnostics* self, size_t instruction,
+                                        bool includeCallers,
                                         bool* instructions) {
   auto* ddd = toDDDiagnostics(self);
   auto* ddsim = ddd->simulationState;
@@ -69,14 +151,37 @@ Result dddiagnosticsGetDataDependencies(Diagnostics* self, size_t instruction,
       instructions, ddsim->interface.getInstructionCount(&ddsim->interface));
   std::set<size_t> toVisit{instruction};
   std::set<size_t> visited;
+
+  // Stores all functions whose callers are unknown (because analysis started
+  // inside them)
+  const std::set<size_t> unknownCallers =
+      includeCallers ? getUnknownCallers(ddsim, instruction)
+                     : std::set<size_t>{};
+
   while (!toVisit.empty()) {
     auto current = *toVisit.begin();
     isDependency[current] = true;
     toVisit.erase(toVisit.begin());
     visited.insert(current);
+
     for (auto dep : ddsim->dataDependencies[current]) {
-      if (visited.find(dep) == visited.end()) {
-        toVisit.insert(dep);
+      const auto depInstruction = dep.first;
+      if (ddsim->instructionTypes[depInstruction] == NOP) {
+        continue; // We don't want variable declarations as dependencies.
+      }
+      if (visited.find(depInstruction) == visited.end()) {
+        toVisit.insert(depInstruction);
+      }
+      if (ddsim->instructionTypes[depInstruction] == CALL) {
+        visitCall(ddsim, depInstruction, dep.second, visited, toVisit);
+      }
+    }
+
+    if (unknownCallers.find(current - 1) != unknownCallers.end()) {
+      for (auto caller : ddsim->functionCallers[current - 1]) {
+        if (visited.find(caller) == visited.end()) {
+          toVisit.insert(caller);
+        }
       }
     }
   }
@@ -95,14 +200,20 @@ Result dddiagnosticsGetInteractions(Diagnostics* self, size_t beforeInstruction,
   while (found) {
     found = false;
     for (auto i = beforeInstruction - 1; i < beforeInstruction; i--) {
+      if (std::find(ddsim->functionDefinitions.begin(),
+                    ddsim->functionDefinitions.end(),
+                    i) != ddsim->functionDefinitions.end()) {
+        break;
+      }
       if (ddsim->instructionTypes[i] != SIMULATE &&
           ddsim->instructionTypes[i] != CALL) {
         continue;
       }
-      auto& targets = ddsim->targetQubits[i];
+
+      auto targets = getTargetVariables(ddsim, i);
       std::set<size_t> targetQubits;
       for (const auto& target : targets) {
-        targetQubits.insert(variableToQubit(ddsim, target));
+        targetQubits.insert(variableToQubitAt(ddsim, target, i).first);
       }
       if (!std::none_of(targetQubits.begin(), targetQubits.end(),
                         [&interactions](size_t elem) {
@@ -150,6 +261,44 @@ size_t dddiagnosticsPotentialErrorCauses(Diagnostics* self, ErrorCause* output,
   return index;
 }
 
+std::set<size_t> getInteractionsAtRuntime(DDDiagnostics* ddd, size_t qubit) {
+  auto* ddsim = ddd->simulationState;
+  std::set<size_t> interactions;
+  interactions.insert(qubit);
+  bool found = true;
+
+  while (found) {
+    found = false;
+
+    for (size_t i = 0; i < ddsim->instructionTypes.size(); i++) {
+      if (ddsim->instructionTypes[i] != SIMULATE) {
+        continue;
+      }
+      if (ddd->actualQubits.find(i) == ddd->actualQubits.end()) {
+        continue;
+      }
+
+      auto& actualQubits = ddd->actualQubits[i];
+      for (const auto& actualQubitVector : actualQubits) {
+        if (!std::none_of(actualQubitVector.begin(), actualQubitVector.end(),
+                          [&interactions](size_t elem) {
+                            return interactions.find(elem) !=
+                                   interactions.end();
+                          })) {
+          for (const auto& target : actualQubitVector) {
+            if (interactions.find(target) == interactions.end()) {
+              found = true;
+            }
+            interactions.insert(target);
+          }
+        }
+      }
+    }
+  }
+
+  return interactions;
+}
+
 size_t tryFindMissingInteraction(DDDiagnostics* diagnostics,
                                  DDSimulationState* state, size_t instruction,
                                  const std::unique_ptr<Assertion>& assertion,
@@ -168,20 +317,18 @@ size_t tryFindMissingInteraction(DDDiagnostics* diagnostics,
                    return variableToQubit(state, target);
                  });
 
-  std::map<size_t, std::vector<uint8_t>> allInteractions;
+  std::map<size_t, std::set<size_t>> allInteractions;
 
   for (size_t i = 0; i < targets.size(); i++) {
-    std::vector<uint8_t> interactions(
-        diagnostics->interface.getNumQubits(&diagnostics->interface));
-    diagnostics->interface.getInteractions(&diagnostics->interface, instruction,
-                                           targetQubits[i],
-                                           toBoolArray(interactions.data()));
-    allInteractions.insert({targetQubits[i], interactions});
+    allInteractions.insert(
+        {targetQubits[i],
+         getInteractionsAtRuntime(diagnostics, targetQubits[i])});
   }
+
   for (size_t i = 0; i < targets.size(); i++) {
     for (size_t j = i + 1; j < targets.size(); j++) {
-      if (allInteractions[targetQubits[i]][targetQubits[j]] == 0 &&
-          allInteractions[targetQubits[j]][targetQubits[i]] == 0) {
+      if (allInteractions[targetQubits[i]].find(targetQubits[j]) ==
+          allInteractions[targetQubits[i]].end()) {
         outputs[index].type = ErrorCauseType::MissingInteraction;
         outputs[index].instruction = instruction;
         index++;
@@ -203,8 +350,9 @@ size_t tryFindZeroControls(DDDiagnostics* diagnostics, size_t instruction,
 
   std::vector<uint8_t> dependencies(
       diagnostics->interface.getInstructionCount(&diagnostics->interface));
-  diagnostics->interface.getDataDependencies(
-      &diagnostics->interface, instruction, toBoolArray(dependencies.data()));
+  diagnostics->interface.getDataDependencies(&diagnostics->interface,
+                                             instruction, true,
+                                             toBoolArray(dependencies.data()));
   auto outputs = Span(output, count);
   size_t index = 0;
 
@@ -270,10 +418,29 @@ Result dddiagnosticsGetZeroControlInstructions(Diagnostics* self,
 void dddiagnosticsOnStepForward(DDDiagnostics* diagnostics,
                                 size_t instruction) {
   auto* ddsim = diagnostics->simulationState;
+  const auto targets = getTargetVariables(ddsim, instruction);
+
+  // Add actual qubits to tracker.
+  if (ddsim->instructionTypes[instruction] == SIMULATE ||
+      ddsim->instructionTypes[instruction] == CALL) {
+    std::vector<size_t> targetQubits(targets.size());
+    std::transform(targets.begin(), targets.end(), targetQubits.begin(),
+                   [&ddsim](const std::string& target) {
+                     return variableToQubit(ddsim, target);
+                   });
+    if (diagnostics->actualQubits.find(instruction) ==
+        diagnostics->actualQubits.end()) {
+      diagnostics->actualQubits[instruction] = std::set<std::vector<size_t>>();
+    }
+    diagnostics->actualQubits[instruction].insert(targetQubits);
+  }
+
+  // Check for zero controls.
   if (ddsim->instructionTypes[instruction] != SIMULATE) {
     return;
   }
-  const auto numQubits = ddsim->interface.getNumQubits(&ddsim->interface);
+  const auto numQubits =
+      diagnostics->interface.getNumQubits(&diagnostics->interface);
   if (numQubits > 16) {
     return;
   }
diff --git a/src/common/parsing/CodePreprocessing.cpp b/src/common/parsing/CodePreprocessing.cpp
index 810ea52..289b70f 100644
--- a/src/common/parsing/CodePreprocessing.cpp
+++ b/src/common/parsing/CodePreprocessing.cpp
@@ -9,14 +9,13 @@
 #include <iterator>
 #include <map>
 #include <memory>
-#include <set>
 #include <string>
 #include <utility>
 #include <vector>
 
 Instruction::Instruction(size_t inputLineNumber, std::string inputCode,
                          std::unique_ptr<Assertion>& inputAssertion,
-                         std::set<std::string> inputTargets, size_t startPos,
+                         std::vector<std::string> inputTargets, size_t startPos,
                          size_t endPos, size_t successor, bool isFuncCall,
                          std::string function, bool inFuncDef, bool isFuncDef,
                          Block inputBlock)
@@ -77,52 +76,82 @@ bool isFunctionDefinition(const std::string& line) {
   return startsWith(trim(line), "gate ");
 }
 
-std::vector<std::string> parseParameters(const std::string& instruction) {
+bool isClassicControlledGate(const std::string& line) {
+  return startsWith(trim(line), "if") &&
+         (line.find('(') != std::string::npos) &&
+         (line.find(')') != std::string::npos);
+}
+
+FunctionDefinition parseFunctionDefinition(const std::string& signature) {
   auto parts = splitString(
-      replaceString(
-          replaceString(replaceString(instruction, ";", " "), "\n", " "), "\t",
-          " "),
-      ' ');
+      replaceString(replaceString(signature, "\n", " "), "\t", " "), ' ');
+  std::string name;
   size_t index = 0;
   for (auto& part : parts) {
     index++;
-    if (!part.empty()) {
+    if (part != "gate" && !part.empty()) {
+      name = part;
       break;
     }
   }
 
   std::string parameterParts;
   for (size_t i = index; i < parts.size(); i++) {
-    if (parts[i].empty()) {
-      continue;
-    }
     parameterParts += parts[i];
   }
   auto parameters = splitString(removeWhitespace(parameterParts), ',');
 
-  return parameters;
+  return {name, parameters};
 }
 
-FunctionDefinition parseFunctionDefinition(const std::string& signature) {
+std::vector<std::string> parseParameters(const std::string& instruction) {
+  if (isFunctionDefinition(instruction)) {
+    const auto fd = parseFunctionDefinition(instruction);
+    return fd.parameters;
+  }
+
+  if (instruction.find("->") != std::string::npos) {
+    // We only add the quantum variable to the measurement's targets.
+    return parseParameters(splitString(instruction, '-')[0]);
+  }
+
+  if (isClassicControlledGate(instruction)) {
+    const auto end = instruction.find(')');
+
+    return parseParameters(
+        instruction.substr(end + 1, instruction.length() - end - 1));
+  }
+
   auto parts = splitString(
-      replaceString(replaceString(signature, "\n", " "), "\t", " "), ' ');
-  std::string name;
+      replaceString(
+          replaceString(replaceString(instruction, ";", " "), "\n", " "), "\t",
+          " "),
+      ' ');
   size_t index = 0;
+  size_t openBrackets = 0;
   for (auto& part : parts) {
     index++;
-    if (part != "gate" && !part.empty()) {
-      name = part;
+    openBrackets +=
+        static_cast<size_t>(std::count(part.begin(), part.end(), '('));
+    openBrackets -=
+        static_cast<size_t>(std::count(part.begin(), part.end(), ')'));
+    if (!part.empty() && openBrackets == 0) {
       break;
     }
   }
 
   std::string parameterParts;
   for (size_t i = index; i < parts.size(); i++) {
+    if (parts[i].empty()) {
+      continue;
+    }
     parameterParts += parts[i];
   }
   auto parameters = splitString(removeWhitespace(parameterParts), ',');
-
-  return {name, parameters};
+  if (parameters.size() == 1 && parameters[0].empty()) {
+    return {};
+  }
+  return parameters;
 }
 
 std::vector<std::string> sweepFunctionNames(const std::string& code) {
@@ -178,9 +207,6 @@ preprocessCode(const std::string& code, size_t startIndex,
     auto tokens = splitString(trimmedLine, ' ');
     auto isAssert = isAssertion(line);
     auto blockPos = line.find("$__block");
-    const auto targetsVector = parseParameters(line);
-    const std::set<std::string> targets(targetsVector.begin(),
-                                        targetsVector.end());
 
     const size_t trueStart = pos + blocksOffset;
 
@@ -199,6 +225,8 @@ preprocessCode(const std::string& code, size_t startIndex,
       line.replace(blockPos, endPos - blockPos + 1, "");
     }
 
+    const auto targets = parseParameters(line);
+
     const size_t trueEnd = end + blocksOffset;
 
     if (isFunctionDefinition(line)) {
@@ -271,20 +299,23 @@ preprocessCode(const std::string& code, size_t startIndex,
   for (auto& instr : instructions) {
     auto vars = parseParameters(instr.code);
     size_t idx = instr.lineNumber - 1;
-    while (!vars.empty() && (instr.lineNumber < instructions.size() ||
-                             idx > instr.lineNumber - instructions.size())) {
-      bool found = false;
+    while (instr.lineNumber != 0 && !vars.empty() &&
+           (instr.lineNumber < instructions.size() ||
+            idx > instr.lineNumber - instructions.size())) {
+      size_t foundIndex = 0;
       for (const auto& var : variableUsages[idx]) {
-        if (std::find(vars.begin(), vars.end(), var) != vars.end()) {
-          found = true;
+        const auto found =
+            std::find_if(vars.begin(), vars.end(), [&var](const auto& v) {
+              return variablesEqual(v, var);
+            });
+        if (found != vars.end()) {
           const auto newEnd = std::remove(vars.begin(), vars.end(), var);
           vars.erase(newEnd, vars.end());
+          instr.dataDependencies.emplace_back(idx, foundIndex);
         }
+        foundIndex++;
       }
-      if (found) {
-        instr.dataDependencies.push_back(idx);
-      }
-      if (idx - 1 == instr.lineNumber - instructions.size()) {
+      if (idx - 1 == instr.lineNumber - instructions.size() || idx == 0) {
         break;
       }
       idx--;
diff --git a/src/common/parsing/Utils.cpp b/src/common/parsing/Utils.cpp
index e6de922..f5f6e3a 100644
--- a/src/common/parsing/Utils.cpp
+++ b/src/common/parsing/Utils.cpp
@@ -55,3 +55,16 @@ std::string removeWhitespace(std::string str) {
   str.erase(std::remove_if(str.begin(), str.end(), ::isspace), str.end());
   return str;
 }
+
+bool variablesEqual(const std::string& v1, const std::string& v2) {
+  if (v1.find('[') != std::string::npos && v2.find('[') != std::string::npos) {
+    return v1 == v2;
+  }
+  if (v1.find('[') != std::string::npos) {
+    return variablesEqual(splitString(v1, '[')[0], v2);
+  }
+  if (v2.find('[') != std::string::npos) {
+    return variablesEqual(splitString(v2, '[')[0], v1);
+  }
+  return v1 == v2;
+}
diff --git a/src/frontend/cli/CliFrontEnd.cpp b/src/frontend/cli/CliFrontEnd.cpp
index e3066a7..0d43053 100644
--- a/src/frontend/cli/CliFrontEnd.cpp
+++ b/src/frontend/cli/CliFrontEnd.cpp
@@ -8,7 +8,6 @@
 #include "backend/diagnostics.h"
 #include "common.h"
 
-#include <array>
 #include <cstddef>
 #include <cstdint>
 #include <iostream>
@@ -127,7 +126,7 @@ void CliFrontEnd::printState(SimulationState* state, size_t inspecting,
     auto* deps = inspectingDependencies.data();
     // NOLINTBEGIN(cppcoreguidelines-pro-type-reinterpret-cast)
     state->getDiagnostics(state)->getDataDependencies(
-        state->getDiagnostics(state), inspecting,
+        state->getDiagnostics(state), inspecting, true,
         reinterpret_cast<bool*>(deps));
     // NOLINTEND(cppcoreguidelines-pro-type-reinterpret-cast)
     uint8_t on = 0;
diff --git a/src/mqt/debug/dap/messages/set_breakpoints_dap_message.py b/src/mqt/debug/dap/messages/set_breakpoints_dap_message.py
index a6633fa..b1a2077 100644
--- a/src/mqt/debug/dap/messages/set_breakpoints_dap_message.py
+++ b/src/mqt/debug/dap/messages/set_breakpoints_dap_message.py
@@ -9,8 +9,6 @@
 if TYPE_CHECKING:
     from .. import DAPServer
 
-# TODO this fails if another file also has breakpoints
-
 
 class SetBreakpointsDAPMessage(DAPMessage):
     """Represents the 'setBreakpoints' DAP request."""
diff --git a/src/mqt/debug/pydebug.pyi b/src/mqt/debug/pydebug.pyi
index 0e9b40d..f6235ea 100644
--- a/src/mqt/debug/pydebug.pyi
+++ b/src/mqt/debug/pydebug.pyi
@@ -89,7 +89,7 @@ class Diagnostics:
     def init(self) -> None: ...
     def get_num_qubits(self) -> int: ...
     def get_instruction_count(self) -> int: ...
-    def get_data_dependencies(self, instruction: int) -> list[int]: ...
+    def get_data_dependencies(self, instruction: int, include_callers: bool = False) -> list[int]: ...
     def get_interactions(self, before_instruction: int, qubit: int) -> list[int]: ...
     def get_zero_control_instructions(self) -> list[int]: ...
     def potential_error_causes(self) -> list[ErrorCause]: ...
diff --git a/src/python/InterfaceBindings.cpp b/src/python/InterfaceBindings.cpp
index a34f8e4..d1a3a2e 100644
--- a/src/python/InterfaceBindings.cpp
+++ b/src/python/InterfaceBindings.cpp
@@ -259,22 +259,24 @@ void bindDiagnostics(py::module& m) {
            [](Diagnostics* self) { return self->getNumQubits(self); })
       .def("get_instruction_count",
            [](Diagnostics* self) { return self->getInstructionCount(self); })
-      .def("get_data_dependencies",
-           [](Diagnostics* self, size_t instruction) {
-             std::vector<uint8_t> instructions(self->getInstructionCount(self));
-             // NOLINTBEGIN(cppcoreguidelines-pro-type-reinterpret-cast)
-             checkOrThrow(self->getDataDependencies(
-                 self, instruction,
-                 reinterpret_cast<bool*>(instructions.data())));
-             // NOLINTEND(cppcoreguidelines-pro-type-reinterpret-cast)
-             std::vector<size_t> result;
-             for (size_t i = 0; i < instructions.size(); i++) {
-               if (instructions[i] != 0) {
-                 result.push_back(i);
-               }
-             }
-             return result;
-           })
+      .def(
+          "get_data_dependencies",
+          [](Diagnostics* self, size_t instruction, bool includeCallers) {
+            std::vector<uint8_t> instructions(self->getInstructionCount(self));
+            // NOLINTBEGIN(cppcoreguidelines-pro-type-reinterpret-cast)
+            checkOrThrow(self->getDataDependencies(
+                self, instruction, includeCallers,
+                reinterpret_cast<bool*>(instructions.data())));
+            // NOLINTEND(cppcoreguidelines-pro-type-reinterpret-cast)
+            std::vector<size_t> result;
+            for (size_t i = 0; i < instructions.size(); i++) {
+              if (instructions[i] != 0) {
+                result.push_back(i);
+              }
+            }
+            return result;
+          },
+          py::arg("instruction"), py::arg("include_callers") = false)
       .def("get_interactions",
            [](Diagnostics* self, size_t beforeInstruction, size_t qubit) {
              std::vector<uint8_t> qubits(self->getNumQubits(self));
diff --git a/test/circuits/diagnose-with-jumps.qasm b/test/circuits/diagnose-with-jumps.qasm
new file mode 100644
index 0000000..8c0d512
--- /dev/null
+++ b/test/circuits/diagnose-with-jumps.qasm
@@ -0,0 +1,26 @@
+gate level_one q0, q1, q2 {                         // 0
+    level_two q0, q1;                               // 1
+    level_two q1, q2;                               // 2
+}                                                   // 3
+
+gate level_two q0, q1 {                             // 4
+    cx q0, q1;                                      // 5
+    level_three_a q0;                               // 6
+
+    level_three_b q1;                               // 7
+}                                                   // 8
+
+gate level_three_a q {                              // 9
+    x q;                                            // 10
+}                                                   // 11
+
+gate level_three_b q {                              // 12
+    z q;                                            // 13
+}                                                   // 14
+
+qreg q[4];                                          // 15
+
+x q[0];                                             // 16
+level_one q[2], q[1], q[0];                         // 17
+
+assert-eq q[0], q[1], q[2], q[3] { qreg x[4]; }     // 18
diff --git a/test/circuits/runtime-interaction.qasm b/test/circuits/runtime-interaction.qasm
new file mode 100644
index 0000000..25efc48
--- /dev/null
+++ b/test/circuits/runtime-interaction.qasm
@@ -0,0 +1,10 @@
+gate gate1 q0, q1 {             // 0
+    x q0;                       // 1
+    x q1;                       // 2
+    assert-ent q0, q1;          // 3
+}                               // 4
+
+qreg q[2];                      // 5
+
+cx q[0], q[1];                  // 6
+gate1 q[0], q[1];               // 7
diff --git a/test/python/resources/diagnosis/data-dependencies-with-callers.qasm b/test/python/resources/diagnosis/data-dependencies-with-callers.qasm
new file mode 100644
index 0000000..03fe2bb
--- /dev/null
+++ b/test/python/resources/diagnosis/data-dependencies-with-callers.qasm
@@ -0,0 +1,15 @@
+qreg q[3];      // 0
+
+gate test q {   // 1
+  x q;          // 2
+}               // 3
+
+x q[0];         // 4
+x q[1];         // 5
+
+test q[0];      // 6
+
+x q[0];         // 7
+x q[2];         // 8
+
+test q[2];      // 9
diff --git a/test/python/test_diagnosis.py b/test/python/test_diagnosis.py
index facba39..b4aa97c 100644
--- a/test/python/test_diagnosis.py
+++ b/test/python/test_diagnosis.py
@@ -44,7 +44,7 @@ def test_data_dependencies_jumps() -> None:
     dependencies = s.get_diagnostics().get_data_dependencies(4)
     assert dependencies == [2, 3, 4]
     dependencies = s.get_diagnostics().get_data_dependencies(9)
-    assert dependencies == [6, 7, 8, 9]
+    assert dependencies == [2, 3, 4, 6, 7, 8, 9]
     destroy_ddsim_simulation_state(s)
 
 
@@ -54,10 +54,13 @@ def test_control_always_zero() -> None:
     s.run_simulation()
     causes = s.get_diagnostics().potential_error_causes()
 
-    assert len(causes) == 1  # once diagnosis can step into jumps, this should be 2
+    assert len(causes) == 2
 
     assert causes[0].type == ErrorCauseType.ControlAlwaysZero
-    assert causes[0].instruction == 12
+    assert causes[0].instruction == 3
+
+    assert causes[1].type == ErrorCauseType.ControlAlwaysZero
+    assert causes[1].instruction == 12
 
 
 def test_missing_interaction() -> None:
@@ -78,3 +81,15 @@ def test_zero_control_listing() -> None:
     s.run_simulation()
     zero_controls = s.get_diagnostics().get_zero_control_instructions()
     assert zero_controls == [3, 12]
+
+
+def test_data_dependencies_with_callers() -> None:
+    """Tests the data dependency analysis with enabling callers."""
+    s = load_instance("data-dependencies-with-callers")
+    s.run_simulation()
+    dependencies = s.get_diagnostics().get_data_dependencies(2, include_callers=True)
+    assert dependencies == [2, 4, 6, 8, 9]
+
+    dependencies = s.get_diagnostics().get_data_dependencies(7, include_callers=True)
+    assert dependencies == [2, 4, 6, 7]
+    # 8 and 9 are not included `test` doesn't have unknown callers in this case, so the analysis won't include all callers.
diff --git a/test/test_custom_code.cpp b/test/test_custom_code.cpp
index 36aa2ff..fbcb732 100644
--- a/test/test_custom_code.cpp
+++ b/test/test_custom_code.cpp
@@ -1,5 +1,6 @@
 #include "backend/dd/DDSimDebug.hpp"
 #include "backend/debug.h"
+#include "backend/diagnostics.h"
 #include "common.h"
 #include "utils_test.hpp"
 
@@ -58,9 +59,9 @@ class CustomCodeTest : public testing::Test {
   }
 };
 
-TEST_F(CustomCodeTest, ClassicControlledOperation) {
+TEST_F(CustomCodeTest, ClassicControlledOperationFalse) {
   loadCode(2, 1,
-           "h q[0];"
+           "z q[0];"
            "cx q[0], q[1];"
            "measure q[0] -> c[0];"
            "if(c==1) x q[1];"
@@ -70,8 +71,24 @@ TEST_F(CustomCodeTest, ClassicControlledOperation) {
   std::array<Complex, 4> amplitudes{};
   Statevector sv{2, 4, amplitudes.data()};
   state->getStateVectorFull(state, &sv);
-  ASSERT_TRUE(complexEquality(amplitudes[0], 1, 0.0) ||
-              complexEquality(amplitudes[1], 1, 0.0));
+  ASSERT_TRUE(complexEquality(amplitudes[0], 1, 0.0));
+
+  ASSERT_EQ(state->stepBackward(state), OK);
+}
+
+TEST_F(CustomCodeTest, ClassicControlledOperationTrue) {
+  loadCode(2, 1,
+           "x q[0];"
+           "cx q[0], q[1];"
+           "measure q[0] -> c[0];"
+           "if(c==1) x q[1];"
+           "if(c==0) z q[1];");
+  ASSERT_EQ(state->runSimulation(state), OK);
+
+  std::array<Complex, 4> amplitudes{};
+  Statevector sv{2, 4, amplitudes.data()};
+  state->getStateVectorFull(state, &sv);
+  ASSERT_TRUE(complexEquality(amplitudes[1], 1, 0.0));
 
   ASSERT_EQ(state->stepBackward(state), OK);
 }
@@ -272,3 +289,124 @@ TEST_F(CustomCodeTest, LargeProgram) {
   ASSERT_EQ(errors, 0);
   ASSERT_EQ(state->getCurrentInstruction(state), 4);
 }
+
+TEST_F(CustomCodeTest, CollectiveGateAsDependency) {
+  loadCode(2, 0, "x q; barrier q[0];");
+  auto* diagnosis = state->getDiagnostics(state);
+  std::array<bool, 4> dependencies{};
+  ASSERT_EQ(
+      diagnosis->getDataDependencies(diagnosis, 3, false, dependencies.data()),
+      OK);
+  ASSERT_EQ(dependencies[0], false);
+  ASSERT_EQ(dependencies[1], false);
+  ASSERT_EQ(dependencies[2], true);
+  ASSERT_EQ(dependencies[3], true);
+}
+
+TEST_F(CustomCodeTest, CollectiveGateAsInteraction) {
+  loadCode(1, 0, "qreg p[1]; cx q, p; assert-ent q[0], p[0];");
+  ASSERT_EQ(state->runSimulation(state), OK);
+  ASSERT_TRUE(state->didAssertionFail(state));
+
+  auto* diagnosis = state->getDiagnostics(state);
+
+  std::array<bool, 2> interactions{};
+  ASSERT_EQ(diagnosis->getInteractions(diagnosis, 4, 0, interactions.data()),
+            OK);
+
+  ASSERT_TRUE(interactions[0]);
+  ASSERT_TRUE(interactions[1]);
+
+  std::array<ErrorCause, 1> causes{};
+  ASSERT_EQ(
+      diagnosis->potentialErrorCauses(diagnosis, causes.data(), causes.size()),
+      1);
+  ASSERT_EQ(causes[0].type, ControlAlwaysZero);
+  ASSERT_EQ(causes[0].instruction, 3);
+}
+
+TEST_F(CustomCodeTest, NonZeroControlsInErrorSearch) {
+  loadCode(2, 0,
+           "gate test q1, q2 { cx q1, q2; } x q[0]; test q[1], q[0]; test "
+           "q[0], q[1]; assert-sup q[0];");
+  auto* diagnosis = state->getDiagnostics(state);
+  ASSERT_EQ(state->runSimulation(state), OK);
+  ASSERT_TRUE(state->didAssertionFail(state));
+  std::array<ErrorCause, 5> errors{};
+  ASSERT_TRUE(diagnosis->potentialErrorCauses(diagnosis, errors.data(),
+                                              errors.size()) == 0);
+}
+
+TEST_F(CustomCodeTest, PaperExampleGrover) {
+  loadCode(3, 3,
+           "gate oracle q0, q1, q2, flag {"
+           "assert-sup q0, q1, q2;"
+           "ccz q1, q2, flag;"
+           "assert-ent q0, q1, q2;"
+           "}"
+           "gate diffusion q0, q1, q2 {"
+           "h q0; h q1; h q2;"
+           "x q0; x q1; x q2;"
+           "ccz q0, q1, q2;"
+           "x q2; x q1; x q0;"
+           "h q2; h q1; h q0;"
+           "}"
+           "qreg flag[1];"
+           "x flag;"
+           "oracle q[0], q[1], q[2], flag;"
+           "diffusion q[0], q[1], q[2];"
+           "assert-eq 0.8, q[0], q[1], q[2] { 0, 0, 0, 0, 0, 0, 0, 1 }"
+           "oracle q[0], q[1], q[2], flag;"
+           "diffusion q[0], q[1], q[2];"
+           "assert-eq 0.9, q[0], q[1], q[2] { 0, 0, 0, 0, 0, 0, 0, 1 }",
+           false, "OPENQASM 2.0;\ninclude \"qelib1.inc\";\n");
+
+  auto* diagnosis = state->getDiagnostics(state);
+  std::array<ErrorCause, 10> causes{};
+
+  ASSERT_EQ(state->runSimulation(state), OK);
+  ASSERT_EQ(state->didAssertionFail(state), true);
+  ASSERT_EQ(state->getCurrentInstruction(state), 5);
+  // We expect no potential errors yet:
+  ASSERT_EQ(
+      diagnosis->potentialErrorCauses(diagnosis, causes.data(), causes.size()),
+      0);
+
+  ASSERT_EQ(state->runSimulation(state), OK);
+  ASSERT_EQ(state->didAssertionFail(state), true);
+  ASSERT_EQ(state->getCurrentInstruction(state), 7);
+  // We expect three potential errors:
+  //   2 missing interactions: q0 <-> q1 and q0 <-> q2
+  //   1 control always zero: q1 & q2 in instruction 6
+  ASSERT_EQ(
+      diagnosis->potentialErrorCauses(diagnosis, causes.data(), causes.size()),
+      3);
+  ASSERT_EQ(causes[0].type, MissingInteraction);
+  ASSERT_EQ(causes[0].instruction, 7);
+  ASSERT_EQ(causes[1].type, MissingInteraction);
+  ASSERT_EQ(causes[1].instruction, 7);
+  ASSERT_EQ(causes[2].type, ControlAlwaysZero);
+  ASSERT_EQ(causes[2].instruction, 6);
+
+  ASSERT_EQ(state->runSimulation(state), OK);
+  ASSERT_EQ(state->didAssertionFail(state), true);
+  ASSERT_EQ(state->getCurrentInstruction(state), 28);
+  // We expect one potential error: Control always zero in instruction 6
+  ASSERT_EQ(
+      diagnosis->potentialErrorCauses(diagnosis, causes.data(), causes.size()),
+      1);
+  ASSERT_EQ(causes[0].type, ControlAlwaysZero);
+  ASSERT_EQ(causes[0].instruction, 6);
+
+  ASSERT_EQ(state->runSimulation(state), OK);
+  ASSERT_EQ(state->didAssertionFail(state), true);
+  ASSERT_EQ(state->getCurrentInstruction(state), 31);
+  // We expect no potential errors, as instruction 6 is no longer always 0
+  ASSERT_EQ(
+      diagnosis->potentialErrorCauses(diagnosis, causes.data(), causes.size()),
+      0);
+
+  ASSERT_EQ(state->runSimulation(state), OK);
+  ASSERT_EQ(state->didAssertionFail(state), false);
+  ASSERT_EQ(state->isFinished(state), true);
+}
diff --git a/test/test_diagnostics.cpp b/test/test_diagnostics.cpp
index d93b430..16f0f58 100644
--- a/test/test_diagnostics.cpp
+++ b/test/test_diagnostics.cpp
@@ -1,6 +1,7 @@
 #include "backend/dd/DDSimDebug.hpp"
 #include "backend/debug.h"
 #include "backend/diagnostics.h"
+#include "common.h"
 #include "utils_test.hpp"
 
 #include <array>
@@ -49,7 +50,8 @@ TEST_F(DiagnosticsTest, DataDependencies) {
     std::vector<uint8_t> dependencies(state->getInstructionCount(state), 0);
     // NOLINTBEGIN(cppcoreguidelines-pro-type-reinterpret-cast)
     diagnostics->getDataDependencies(
-        diagnostics, instruction, reinterpret_cast<bool*>(dependencies.data()));
+        diagnostics, instruction, false,
+        reinterpret_cast<bool*>(dependencies.data()));
     // NOLINTEND(cppcoreguidelines-pro-type-reinterpret-cast)
     std::set<size_t> dependenciesSet;
     for (size_t i = 0; i < dependencies.size(); ++i) {
@@ -171,3 +173,98 @@ TEST_F(DiagnosticsTest, ZeroControlsWithJumps) {
     ASSERT_FALSE(zeroControls.at(i) ^ (i == 3 || i == 12));
   }
 }
+
+TEST_F(DiagnosticsTest, DataDependenciesWithJumps) {
+  loadFromFile("diagnose-with-jumps");
+  const std::map<size_t, std::set<size_t>> expected = {
+      {1, {1}},
+      {2, {2, 13, 7, 5, 1}},
+      {3, {3}},
+
+      {5, {5}},
+      {6, {6, 5}},
+      {7, {7, 5}},
+      {8, {8}},
+
+      {10, {10}},
+      {13, {13}},
+      {11, {11}},
+      {9, {9}},
+      {14, {14}},
+      {12, {12}},
+
+      {15, {15}},
+
+      {16, {16}},
+
+      {17, {17, 16}},
+
+      {18, {18, 13, 10, 7, 6, 5, 2, 1, 17, 16}}};
+
+  for (const auto& pair : expected) {
+    std::vector<uint8_t> dependencies(state->getInstructionCount(state), 0);
+    // NOLINTBEGIN(cppcoreguidelines-pro-type-reinterpret-cast)
+    diagnostics->getDataDependencies(
+        diagnostics, pair.first, false,
+        reinterpret_cast<bool*>(dependencies.data()));
+    // NOLINTEND(cppcoreguidelines-pro-type-reinterpret-cast)
+    std::set<size_t> dependenciesSet;
+    for (size_t i = 0; i < dependencies.size(); ++i) {
+      if (dependencies[i] != 0) {
+        dependenciesSet.insert(i);
+      }
+    }
+    ASSERT_EQ(dependenciesSet, pair.second)
+        << "Failed for instruction " << pair.first;
+  }
+}
+
+TEST_F(DiagnosticsTest, InteractionsWithJumps) {
+  loadFromFile("diagnose-with-jumps");
+
+  const std::map<std::pair<size_t, size_t>, std::set<size_t>> expected = {
+      {{1, 0}, {0}},        {{1, 1}, {1}},        {{1, 2}, {2}},
+      {{2, 0}, {0, 1}},     {{2, 1}, {0, 1}},     {{2, 2}, {2}},
+
+      {{5, 0}, {0}},        {{6, 0}, {1, 0}},     {{7, 1}, {0, 1}},
+
+      {{10, 0}, {0}},
+
+      {{17, 0}, {0}},       {{18, 0}, {1, 2, 0}}, {{18, 1}, {0, 2, 1}},
+      {{18, 2}, {0, 1, 2}}, {{18, 3}, {3}}};
+
+  for (const auto& pair : expected) {
+    std::vector<uint8_t> interactions(state->getNumQubits(state), 0);
+    // NOLINTBEGIN(cppcoreguidelines-pro-type-reinterpret-cast)
+    diagnostics->getInteractions(diagnostics, pair.first.first,
+                                 pair.first.second,
+                                 reinterpret_cast<bool*>(interactions.data()));
+    // NOLINTEND(cppcoreguidelines-pro-type-reinterpret-cast)
+    std::set<size_t> interactionsSet;
+    for (size_t i = 0; i < interactions.size(); ++i) {
+      if (interactions[i] != 0) {
+        interactionsSet.insert(i);
+      }
+    }
+    ASSERT_EQ(interactionsSet, pair.second)
+        << "Failed for instruction " << pair.first.first << " qubit "
+        << pair.first.second;
+  }
+}
+
+TEST_F(DiagnosticsTest, RuntimeInteractions) {
+  loadFromFile("runtime-interaction");
+  ASSERT_EQ(state->runSimulation(state), OK);
+  ASSERT_TRUE(state->didAssertionFail(state));
+  ASSERT_EQ(state->getCurrentInstruction(state), 3);
+
+  std::array<ErrorCause, 10> errors{};
+  ASSERT_EQ(diagnostics->potentialErrorCauses(diagnostics, errors.data(), 10),
+            1);
+
+  // Interaction happens outside of the instruction, so we don't expect a
+  // missing interaction.
+
+  ASSERT_EQ(errors[0].type, ErrorCauseType::ControlAlwaysZero);
+  ASSERT_EQ(errors[0].instruction, 6);
+}