diff --git a/deerlab/model.py b/deerlab/model.py
index a8f2cca3..81f5fac2 100644
--- a/deerlab/model.py
+++ b/deerlab/model.py
@@ -1018,7 +1018,7 @@ def decorator(func):
 #==============================================================================================
 @insert_snlls_optionals_docstrings()
 def fit(model_, y, *constants, par0=None, penalties=None, bootstrap=0, noiselvl=None, mask=None, weights=None,
-                regparam='aic',reg='auto',regparamrange=None,**kwargs):
+                regparam='aic',reg='auto',regparamrange=None, bootcores=1,**kwargs):
     r"""
     Fit the model(s) to the dataset(s)
 
@@ -1051,6 +1051,10 @@ def fit(model_, y, *constants, par0=None, penalties=None, bootstrap=0, noiselvl=
     bootstrap : scalar, optional,
         Bootstrap samples for uncertainty quantification. If ``bootstrap>0``, the uncertainty quantification will be 
         performed via the boostrapping method with based on the number of samples specified as the argument.
+    
+    bootcores : scalar, optional
+        Number of CPU cores/processes for parallelization of the bootstrap uncertainty quantification. If ``cores=1`` no parallel 
+        computing is used. If ``cores=-1`` all available CPUs are used. The default is one core (no parallelization).
 
     reg : boolean or string, optional
         Determines the use of regularization on the solution of the linear problem.
@@ -1264,7 +1268,7 @@ def bootstrap_fcn(ysim):
             if not isinstance(fit.model,list): fit.model = [fit.model]
             return (fit.param,*fit.model)
         # Bootstrapped uncertainty quantification
-        param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap,noiselvl=noiselvl)
+        param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap,noiselvl=noiselvl,cores=bootcores)
         # Include information on the boundaries for better uncertainty estimates
         paramlb = model._vecsort(model._getvector('lb'))[np.concatenate(param_idx)] 
         paramub = model._vecsort(model._getvector('ub'))[np.concatenate(param_idx)]