diff --git a/Cargo.toml b/Cargo.toml
index 92965d9a..fd48b7f1 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -7,27 +7,51 @@ description = "Rust library for solving analytic and ode-defined pharmacometric
 license = "GPL-3.0"
 documentation = "https://lapkb.github.io/pharmsol/"
 
+
 [features]
-default = []
+
+
+default = [
+    "native"
+]
+
+# `native` now contains only crates that are known to break or behave poorly
+# on wasm (threads, native-optimized code, or platform-specific backends).
+native = [
+    "dep:rayon",
+    "ndarray/rayon",
+    "dep:argmin",
+    "dep:argmin-math",
+    "dep:getrandom",
+    "dep:statrs",
+]
+
+wasm = []
+
 exa = ["libloading"]
 
 [dependencies]
 cached = { version = "0.56.0" }
 csv = "1.3.0"
-diffsol = "=0.7.0"
+diffsol = { version = "=0.7.0",  features = ["faer", "nalgebra"] }
 libloading = { version = "0.8.6", optional = true, features = [] }
-nalgebra = "0.34.1"
-ndarray = { version = "0.16.1", features = ["rayon"] }
-rand = "0.9.0"
-rand_distr = "0.5.0"
-rayon = "1.10.0"
+nalgebra = { version = "0.34.1" }
+ndarray = { version = "0.16.1"}
+rand = { version = "0.9.0"}
+rand_distr = { version = "0.5.0"}
+rayon = { version = "1.10.0", optional = true }
 serde = { version = "1.0.201", features = ["derive"] }
 serde_json = "1.0.117"
-statrs = "0.18.0"
+statrs = { version = "0.18.0", optional = true }
 thiserror = "2.0.11"
-argmin = "0.11.0"
-argmin-math = "0.5.1"
+argmin = { version = "0.11.0", optional = true }
+argmin-math = { version = "0.5.1", optional = true }
 tracing = "0.1.41"
+getrandom = { version = "0.3", optional = true }
+
+[target.'cfg(target_arch = "wasm32")'.dependencies]
+
+getrandom = { version = "0.3", features = ["wasm_js"] }
 
 [dev-dependencies]
 criterion = { version = "0.7.0", features = ["html_reports"] }
diff --git a/src/lib.rs b/src/lib.rs
index 58f01d8b..f9522cac 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -10,7 +10,10 @@ pub use crate::data::builder::SubjectBuilderExt;
 pub use crate::data::Interpolation::*;
 pub use crate::data::*;
 pub use crate::equation::*;
+#[cfg(feature = "native")]
 pub use crate::optimize::effect::get_e2;
+
+#[cfg(feature = "native")]
 pub use crate::optimize::spp::SppOptimizer;
 pub use crate::simulator::equation::{self, ODE};
 pub use error::PharmsolError;
diff --git a/src/optimize/effect.rs b/src/optimize/effect.rs
index f132393d..256babd5 100644
--- a/src/optimize/effect.rs
+++ b/src/optimize/effect.rs
@@ -1,8 +1,10 @@
+#[cfg(feature = "native")]
 use argmin::{
     core::{CostFunction, Executor, TerminationReason, TerminationStatus},
     solver::neldermead::NelderMead,
 };
 
+#[cfg(feature = "native")]
 #[derive(Debug, Clone)]
 struct BestM0 {
     a: f64,
@@ -14,6 +16,7 @@ struct BestM0 {
 }
 
 /// We'll optimize over y = ln(xm), so Param = f64 (the log of xm)
+#[cfg(feature = "native")]
 impl CostFunction for BestM0 {
     type Param = f64; // this is ln(xm)
     type Output = f64;
@@ -59,12 +62,14 @@ impl CostFunction for BestM0 {
     }
 }
 
+#[cfg(feature = "native")]
 #[derive(thiserror::Error, Debug)]
 pub enum BestM0Error {
     #[error("Optimization error: {0}")]
     OptimizationError(String),
 }
 
+#[cfg(feature = "native")]
 impl BestM0 {
     /// start and step are in log-space (ln(x))
     fn get_best(&self, start_log: f64, step_log: f64) -> Result<(f64, f64, bool), BestM0Error> {
@@ -111,8 +116,7 @@ impl BestM0 {
     }
 }
 
-/// find_m0 left largely as-is, but consider returning Result<f64>
-/// Keep in mind it expects a,b,h1,h2 in valid ranges
+#[cfg(feature = "native")]
 fn find_m0(afinal: f64, b: f64, alpha: f64, h1: f64, h2: f64) -> f64 {
     let noint = 1000;
     let del_a = afinal / (noint as f64);
@@ -150,6 +154,7 @@ fn find_m0(afinal: f64, b: f64, alpha: f64, h1: f64, h2: f64) -> f64 {
     xm
 }
 
+#[cfg(feature = "native")]
 pub fn get_e2(a: f64, b: f64, w: f64, h1: f64, h2: f64, alpha_s: f64) -> f64 {
     // trivial cases
     if a.abs() < 1.0e-12 && b.abs() < 1.0e-12 {
@@ -237,3 +242,11 @@ pub fn get_e2(a: f64, b: f64, w: f64, h1: f64, h2: f64, alpha_s: f64) -> f64 {
         }
     }
 }
+
+// Non-native (wasm) stub implementation
+#[cfg(not(feature = "native"))]
+pub fn get_e2(_a: f64, _b: f64, _w: f64, _h1: f64, _h2: f64, _alpha_s: f64) -> f64 {
+    // Provide a conservative default for wasm builds. This is a best-effort
+    // stub — accurate behavior requires the native optimizer.
+    0.0
+}
diff --git a/src/optimize/spp.rs b/src/optimize/spp.rs
index ac19bd3b..ec0b365c 100644
--- a/src/optimize/spp.rs
+++ b/src/optimize/spp.rs
@@ -1,12 +1,16 @@
+#[cfg(feature = "native")]
 use argmin::{
     core::{CostFunction, Error, Executor},
     solver::neldermead::NelderMead,
 };
 
+#[cfg(feature = "native")]
 use ndarray::{Array1, Axis};
 
+#[cfg(feature = "native")]
 use crate::{prelude::simulator::psi, Data, Equation, ErrorModels};
 
+#[cfg(feature = "native")]
 pub struct SppOptimizer<'a, E: Equation> {
     equation: &'a E,
     data: &'a Data,
@@ -14,6 +18,7 @@ pub struct SppOptimizer<'a, E: Equation> {
     pyl: &'a Array1<f64>,
 }
 
+#[cfg(feature = "native")]
 impl<E: Equation> CostFunction for SppOptimizer<'_, E> {
     type Param = Vec<f64>;
     type Output = f64;
@@ -41,6 +46,7 @@ impl<E: Equation> CostFunction for SppOptimizer<'_, E> {
     }
 }
 
+#[cfg(feature = "native")]
 impl<'a, E: Equation> SppOptimizer<'a, E> {
     pub fn new(
         equation: &'a E,
@@ -66,6 +72,7 @@ impl<'a, E: Equation> SppOptimizer<'a, E> {
     }
 }
 
+#[cfg(feature = "native")]
 fn create_initial_simplex(initial_point: &[f64]) -> Vec<Vec<f64>> {
     let num_dimensions = initial_point.len();
     let perturbation_percentage = 0.008;
@@ -91,3 +98,39 @@ fn create_initial_simplex(initial_point: &[f64]) -> Vec<Vec<f64>> {
 
     vertices
 }
+
+// WASM / non-native stubs
+#[cfg(not(feature = "native"))]
+use ndarray::Array1;
+
+#[cfg(not(feature = "native"))]
+use crate::{Data, Equation, ErrorModels};
+
+#[cfg(not(feature = "native"))]
+pub struct SppOptimizer<'a, E> {
+    _phantom: std::marker::PhantomData<&'a E>,
+}
+
+#[cfg(not(feature = "native"))]
+impl<'a, E: Equation> SppOptimizer<'a, E> {
+    pub fn new(
+        _equation: &'a E,
+        _data: &'a Data,
+        _sig: &'a ErrorModels,
+        _pyl: &'a Array1<f64>,
+    ) -> Self {
+        SppOptimizer {
+            _phantom: std::marker::PhantomData,
+        }
+    }
+
+    pub fn optimize_point(
+        self,
+        _spp: Array1<f64>,
+    ) -> Result<Array1<f64>, Box<dyn std::error::Error>> {
+        Err(Box::new(std::io::Error::new(
+            std::io::ErrorKind::Other,
+            "SppOptimizer is disabled when feature `native` is not enabled",
+        )))
+    }
+}
diff --git a/src/simulator/equation/sde/mod.rs b/src/simulator/equation/sde/mod.rs
index 85168378..c2268a3f 100644
--- a/src/simulator/equation/sde/mod.rs
+++ b/src/simulator/equation/sde/mod.rs
@@ -4,6 +4,8 @@ use diffsol::{NalgebraContext, Vector};
 use nalgebra::DVector;
 use ndarray::{concatenate, Array2, Axis};
 use rand::{rng, Rng};
+
+#[cfg(feature = "native")]
 use rayon::prelude::*;
 
 use cached::proc_macro::cached;
@@ -142,9 +144,18 @@ impl State for Vec<DVector<f64>> {
     /// * `input` - Index of the input compartment
     /// * `amount` - Amount to add to the compartment
     fn add_bolus(&mut self, input: usize, amount: f64) {
-        self.par_iter_mut().for_each(|particle| {
-            particle[input] += amount;
-        });
+        #[cfg(feature = "native")]
+        {
+            self.par_iter_mut().for_each(|particle| {
+                particle[input] += amount;
+            });
+        }
+        #[cfg(not(feature = "native"))]
+        {
+            self.iter_mut().for_each(|particle| {
+                particle[input] += amount;
+            });
+        }
     }
 }
 
@@ -239,20 +250,40 @@ impl EquationPriv for SDE {
         ti: f64,
         tf: f64,
     ) -> Result<(), PharmsolError> {
-        state.par_iter_mut().for_each(|particle| {
-            *particle = simulate_sde_event(
-                &self.drift,
-                &self.diffusion,
-                particle.clone().into(),
-                support_point,
-                covariates,
-                infusions,
-                ti,
-                tf,
-            )
-            .inner()
-            .clone();
-        });
+        #[cfg(feature = "native")]
+        {
+            state.par_iter_mut().for_each(|particle| {
+                *particle = simulate_sde_event(
+                    &self.drift,
+                    &self.diffusion,
+                    particle.clone().into(),
+                    support_point,
+                    covariates,
+                    infusions,
+                    ti,
+                    tf,
+                )
+                .inner()
+                .clone();
+            });
+        }
+        #[cfg(not(feature = "native"))]
+        {
+            for particle in state.iter_mut() {
+                *particle = simulate_sde_event(
+                    &self.drift,
+                    &self.diffusion,
+                    particle.clone().into(),
+                    support_point,
+                    covariates,
+                    infusions,
+                    ti,
+                    tf,
+                )
+                .inner()
+                .clone();
+            }
+        }
         Ok(())
     }
     fn nparticles(&self) -> usize {
@@ -275,17 +306,34 @@ impl EquationPriv for SDE {
         output: &mut Self::P,
     ) -> Result<(), PharmsolError> {
         let mut pred = vec![Prediction::default(); self.nparticles];
-        pred.par_iter_mut().enumerate().for_each(|(i, p)| {
-            let mut y = V::zeros(self.get_nouteqs(), NalgebraContext);
-            (self.out)(
-                &x[i].clone().into(),
-                &V::from_vec(support_point.clone(), NalgebraContext),
-                observation.time(),
-                covariates,
-                &mut y,
-            );
-            *p = observation.to_prediction(y[observation.outeq()], x[i].as_slice().to_vec());
-        });
+        #[cfg(feature = "native")]
+        {
+            pred.par_iter_mut().enumerate().for_each(|(i, p)| {
+                let mut y = V::zeros(self.get_nouteqs(), NalgebraContext);
+                (self.out)(
+                    &x[i].clone().into(),
+                    &V::from_vec(support_point.clone(), NalgebraContext),
+                    observation.time(),
+                    covariates,
+                    &mut y,
+                );
+                *p = observation.to_prediction(y[observation.outeq()], x[i].as_slice().to_vec());
+            });
+        }
+        #[cfg(not(feature = "native"))]
+        {
+            for (i, p) in pred.iter_mut().enumerate() {
+                let mut y = V::zeros(self.get_nouteqs(), NalgebraContext);
+                (self.out)(
+                    &x[i].clone().into(),
+                    &V::from_vec(support_point.clone(), NalgebraContext),
+                    observation.time(),
+                    covariates,
+                    &mut y,
+                );
+                *p = observation.to_prediction(y[observation.outeq()], x[i].as_slice().to_vec());
+            }
+        }
         let out = Array2::from_shape_vec((self.nparticles, 1), pred.clone())?;
         *output = concatenate(Axis(1), &[output.view(), out.view()]).unwrap();
         //e = y[t] .- x[:,1]
diff --git a/src/simulator/likelihood/mod.rs b/src/simulator/likelihood/mod.rs
index b4401775..1c9b049c 100644
--- a/src/simulator/likelihood/mod.rs
+++ b/src/simulator/likelihood/mod.rs
@@ -6,8 +6,14 @@ use crate::{
     Predictions,
 };
 use ndarray::{Array2, Axis, ShapeBuilder};
+
+#[cfg(feature = "native")]
 use rayon::prelude::*;
+
+#[cfg(feature = "native")]
 use statrs::distribution::ContinuousCDF;
+
+#[cfg(feature = "native")]
 use statrs::distribution::Normal;
 
 mod progress;
@@ -108,12 +114,40 @@ impl SubjectPredictions {
 fn normpdf(obs: f64, pred: f64, sigma: f64) -> f64 {
     (FRAC_1_SQRT_2PI / sigma) * (-((obs - pred) * (obs - pred)) / (2.0 * sigma * sigma)).exp()
 }
+
+#[cfg(feature = "native")]
 #[inline(always)]
 fn normcdf(obs: f64, pred: f64, sigma: f64) -> Result<f64, ErrorModelError> {
     let norm = Normal::new(pred, sigma).map_err(|_| ErrorModelError::NegativeSigma)?;
     Ok(norm.cdf(obs))
 }
 
+// wasm-safe approximation of the normal CDF (avoids `statrs`)
+#[cfg(not(feature = "native"))]
+#[inline(always)]
+fn normcdf(obs: f64, pred: f64, sigma: f64) -> Result<f64, ErrorModelError> {
+    if sigma <= 0.0 {
+        return Err(ErrorModelError::NegativeSigma);
+    }
+    // Use a numeric approximation for erf (Abramowitz-Stegun)
+    fn erf(x: f64) -> f64 {
+        let a1 = 0.254829592;
+        let a2 = -0.284496736;
+        let a3 = 1.421413741;
+        let a4 = -1.453152027;
+        let a5 = 1.061405429;
+        let p = 0.3275911;
+        let sign = if x < 0.0 { -1.0 } else { 1.0 };
+        let x = x.abs();
+        let t = 1.0 / (1.0 + p * x);
+        let y = 1.0 - (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t) * (-x * x).exp();
+        sign * y
+    }
+
+    let z = (obs - pred) / (sigma * std::f64::consts::SQRT_2);
+    Ok(0.5 * (1.0 + erf(z)))
+}
+
 impl From<Vec<Prediction>> for SubjectPredictions {
     fn from(predictions: Vec<Prediction>) -> Self {
         Self {
@@ -182,15 +216,42 @@ pub fn psi(
         None
     };
 
-    let result: Result<(), PharmsolError> = psi
-        .axis_iter_mut(Axis(0))
-        .into_par_iter()
-        .enumerate()
-        .try_for_each(|(i, mut row)| {
-            row.axis_iter_mut(Axis(0))
+    let result: Result<(), PharmsolError> = {
+        #[cfg(feature = "native")]
+        {
+            psi.axis_iter_mut(Axis(0))
                 .into_par_iter()
                 .enumerate()
-                .try_for_each(|(j, mut element)| {
+                .try_for_each(|(i, mut row)| {
+                    row.axis_iter_mut(Axis(0))
+                        .into_par_iter()
+                        .enumerate()
+                        .try_for_each(|(j, mut element)| {
+                            let subject = subjects.get(i).unwrap();
+                            match equation.estimate_likelihood(
+                                subject,
+                                support_points.row(j).to_vec().as_ref(),
+                                error_models,
+                                cache,
+                            ) {
+                                Ok(likelihood) => {
+                                    element.fill(likelihood);
+                                    if let Some(ref tracker) = progress_tracker {
+                                        tracker.inc();
+                                    }
+                                }
+                                Err(e) => return Err(e),
+                            };
+                            Ok(())
+                        })
+                })
+        }
+
+        #[cfg(not(feature = "native"))]
+        {
+            let mut outer_result: Result<(), PharmsolError> = Ok(());
+            for (i, mut row) in psi.axis_iter_mut(Axis(0)).enumerate() {
+                for (j, mut element) in row.axis_iter_mut(Axis(0)).enumerate() {
                     let subject = subjects.get(i).unwrap();
                     match equation.estimate_likelihood(
                         subject,
@@ -204,11 +265,19 @@ pub fn psi(
                                 tracker.inc();
                             }
                         }
-                        Err(e) => return Err(e),
+                        Err(e) => {
+                            outer_result = Err(e);
+                            break;
+                        }
                     };
-                    Ok(())
-                })
-        });
+                }
+                if outer_result.is_err() {
+                    break;
+                }
+            }
+            outer_result
+        }
+    };
 
     if let Some(tracker) = progress_tracker {
         tracker.finish();