diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root
new file mode 120000
index 000000000..945c9b46d
--- /dev/null
+++ b/_codeql_detected_source_root
@@ -0,0 +1 @@
+.
\ No newline at end of file
diff --git a/src/expm.cpp b/src/expm.cpp
index e75791822..970163a56 100644
--- a/src/expm.cpp
+++ b/src/expm.cpp
@@ -310,12 +310,13 @@ int meOnly(int cSub, double *yc_, double *yp_, double tp, double tf, double tcov
 //'        inductive linearization
 //' @name rxIndLin_
 //' @noRd
-extern "C" int indLin(int cSub, rx_solving_options *op, double tp, double *yp_, double tf,
+extern "C" int indLin(int cSub, rx_solving_options *op, rx_solving_options_ind *ind, double tp, double *yp_, double tf,
 		      double *InfusionRate_, int *on_,
 		      t_ME ME, t_IndF  IndF){
   int neq = op->neq;
-  double *rtol=op->rtol2;
-  double *atol=op->atol2;
+  // Use per-individual tolerance arrays if available
+  double *rtol = (ind != NULL && ind->rtol2 != NULL) ? ind->rtol2 : op->rtol2;
+  double *atol = (ind != NULL && ind->atol2 != NULL) ? ind->atol2 : op->atol2;
   int maxsteps=op->mxstep;
   int doIndLin=op->doIndLin;
   // int indLinPerterb=10;
diff --git a/src/par_solve.cpp b/src/par_solve.cpp
index b23e6b052..c77993658 100644
--- a/src/par_solve.cpp
+++ b/src/par_solve.cpp
@@ -732,7 +732,7 @@ static const char *err_msg_ls[] =
 //dummy solout fn
 extern "C" void solout(long int nr, double t_old, double t, double *y, int *nptr, int *irtrn){}
 
-extern "C" int indLin(int cSub, rx_solving_options *op, double tp, double *yp_, double tf,
+extern "C" int indLin(int cSub, rx_solving_options *op, rx_solving_options_ind *ind, double tp, double *yp_, double tf,
                       double *InfusionRate_, int *on_,
                       t_ME ME, t_IndF  IndF);
 
@@ -760,7 +760,7 @@ static inline void solveWith1Pt(int *neq,
     case 3:
       if (!isSameTime(xout, xp)) {
         preSolve(op, ind, xp, xout, yp);
-        idid = indLin(ind->id, op, xp, yp, xout, ind->InfusionRate, ind->on,
+        idid = indLin(ind->id, op, ind, xp, yp, xout, ind->InfusionRate, ind->on,
                       ME, IndF);
       }
       if (idid <= 0) {
@@ -2251,7 +2251,7 @@ extern "C" void ind_indLin0(rx_solve *rx, rx_solving_options *op, int solveid,
         badSolveExit(i);
       } else {
         preSolve(op, ind, xoutp, xout, yp);
-        idid = indLin(solveid, op, xoutp, yp, xout, ind->InfusionRate, ind->on,
+        idid = indLin(solveid, op, ind, xoutp, yp, xout, ind->InfusionRate, ind->on,
                       ME, IndF);
         xoutp=xout;
         postSolve(neq, &idid, rc, &i, yp, NULL, 0, true, ind, op, rx);
@@ -2368,6 +2368,13 @@ extern "C" void ind_liblsoda0(rx_solve *rx, rx_solving_options *op, struct lsoda
     ctx = NULL;
     return;
   }
+  
+  // Use per-individual tolerance arrays if available
+  if (ind->rtol2 != NULL && ind->atol2 != NULL) {
+    opt.rtol = ind->rtol2;
+    opt.atol = ind->atol2;
+  }
+  
   nx = ind->n_all_times;
   BadDose = ind->BadDose;
   InfusionRate = ind->InfusionRate;
diff --git a/src/rxData.cpp b/src/rxData.cpp
index 2b9aabc88..600a7538c 100644
--- a/src/rxData.cpp
+++ b/src/rxData.cpp
@@ -1540,6 +1540,12 @@ struct rx_globals {
   // time per thread
   double *timeThread = NULL;
 
+  // Per-thread tolerance arrays for thread-safe atolRtolFactor_
+  double *gatol2Thread = NULL;
+  double *grtol2Thread = NULL;
+  double *gssAtolThread = NULL;
+  double *gssRtolThread = NULL;
+
   bool alloc=false;
 };
 
@@ -1633,6 +1639,12 @@ extern "C" void _setIndPointersByThread(rx_solving_options_ind *ind) {
   ind->timeThread = _globals.timeThread + rx->maxAllTimes*omp_get_thread_num();
   ind->llikSave = _globals.gLlikSave + op->nLlik*rxLlikSaveSize*omp_get_thread_num();
   ind->lhs = _globals.glhs+op->nlhs*omp_get_thread_num();
+  
+  // Assign per-thread tolerance arrays
+  ind->atol2 = _globals.gatol2Thread + op->neq*omp_get_thread_num();
+  ind->rtol2 = _globals.grtol2Thread + op->neq*omp_get_thread_num();
+  ind->ssAtol = _globals.gssAtolThread + op->neq*omp_get_thread_num();
+  ind->ssRtol = _globals.gssRtolThread + op->neq*omp_get_thread_num();
 }
 
 extern "C" void setZeroMatrix(int which) {
@@ -1655,12 +1667,23 @@ double maxAtolRtolFactor = 0.1;
 void atolRtolFactor_(double factor){
   rx_solve* rx = getRxSolve_();
   rx_solving_options* op = rx->op;
-  for (int i = op->neq;i--;){
-    _globals.grtol2[i] = min2(_globals.grtol2[i]*factor, maxAtolRtolFactor);
-    _globals.gatol2[i] = min2(_globals.gatol2[i]*factor, maxAtolRtolFactor);
-  }
-  op->ATOL = min2(op->ATOL*factor, maxAtolRtolFactor);
-  op->RTOL = min2(op->RTOL*factor, maxAtolRtolFactor);
+  
+  // Get thread-specific tolerance arrays
+  int threadId = omp_get_thread_num();
+  double *atol2 = _globals.gatol2Thread + op->neq*threadId;
+  double *rtol2 = _globals.grtol2Thread + op->neq*threadId;
+  double *ssAtol = _globals.gssAtolThread + op->neq*threadId;
+  double *ssRtol = _globals.gssRtolThread + op->neq*threadId;
+  
+  // Modify thread-specific tolerances
+  for (int i = op->neq; i--;){
+    atol2[i] = min2(atol2[i]*factor, maxAtolRtolFactor);
+    rtol2[i] = min2(rtol2[i]*factor, maxAtolRtolFactor);
+    ssAtol[i] = min2(ssAtol[i]*factor, maxAtolRtolFactor);
+    ssRtol[i] = min2(ssRtol[i]*factor, maxAtolRtolFactor);
+  }
+  // Note: op->ATOL and op->RTOL are not modified to avoid race conditions
+  // The per-thread tolerance arrays are used by the solvers
 }
 
 extern "C" double * getAol(int n, double atol){
@@ -5256,10 +5279,11 @@ SEXP rxSolve_(const RObject &obj, const List &rxControl,
     int n8 = rx->maxAllTimes*op->cores;
     int n9 = (op->numLinSens+op->numLin)*op->cores;
     int n10 = (op->neq)*op->cores;
+    int n11 = 4*op->neq*op->cores; // Per-thread tolerance arrays (atol2, rtol2, ssAtol, ssRtol)
     int nlin = (rx->linB)* 7* rx->nsub * rx->nsim;
     if (_globals.gsolve != NULL) free(_globals.gsolve);
     _globals.gsolve = (double*)calloc(nlin+n0+3*nsave+n2+ n4+n5_c+n6+ n7 + n8 +
-                                      n9 + n10 +
+                                      n9 + n10 + n11 +
                                       5*op->neq + 4*n3a_c + nllik_c,
                                       sizeof(double));// [n0]
 #ifdef rxSolveT
@@ -5298,7 +5322,11 @@ SEXP rxSolve_(const RObject &obj, const List &rxControl,
     _globals.gIndSim   = _globals.gCurDoseS + n3a_c;// [n7]
     _globals.gLinSave  = _globals.gIndSim + n7; // [n9]
     _globals.gLinDummy = _globals.gLinSave + n9; // [n10]
-    _globals.timeThread = _globals.gLinDummy + n10;
+    _globals.timeThread = _globals.gLinDummy + n10; // [n8]
+    _globals.gatol2Thread = _globals.timeThread + n8; // [op->neq * op->cores]
+    _globals.grtol2Thread = _globals.gatol2Thread + op->neq*op->cores; // [op->neq * op->cores]
+    _globals.gssAtolThread = _globals.grtol2Thread + op->neq*op->cores; // [op->neq * op->cores]
+    _globals.gssRtolThread = _globals.gssAtolThread + op->neq*op->cores; // [op->neq * op->cores]
     std::fill_n(rx->ypNA, op->neq, NA_REAL);
 
     std::fill_n(&_globals.gatol2[0],op->neq, atolNV[0]);
@@ -5310,6 +5338,19 @@ SEXP rxSolve_(const RObject &obj, const List &rxControl,
     std::fill_n(&_globals.gssRtol[0],op->neq, rtolNVss[0]);
     std::copy(atolNVss.begin(), atolNVss.begin() + min2(op->neq, atolNVss.size()), &_globals.gssAtol[0]);
     std::copy(rtolNVss.begin(), rtolNVss.begin() + min2(op->neq, rtolNVss.size()), &_globals.gssRtol[0]);
+    
+    // Initialize per-thread tolerance arrays from global values
+    for (int core = 0; core < op->cores; core++) {
+      std::copy(&_globals.gatol2[0], &_globals.gatol2[0] + op->neq, 
+                _globals.gatol2Thread + core*op->neq);
+      std::copy(&_globals.grtol2[0], &_globals.grtol2[0] + op->neq,
+                _globals.grtol2Thread + core*op->neq);
+      std::copy(&_globals.gssAtol[0], &_globals.gssAtol[0] + op->neq,
+                _globals.gssAtolThread + core*op->neq);
+      std::copy(&_globals.gssRtol[0], &_globals.gssRtol[0] + op->neq,
+                _globals.gssRtolThread + core*op->neq);
+    }
+    
     op->atol2 = &_globals.gatol2[0];
     op->rtol2 = &_globals.grtol2[0];
     op->ssAtol = _globals.gssAtol;