sc_neurocore_engine/
dna.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later
2// Commercial license available
3// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
4// © Code 2020–2026 Miroslav Šotek. All rights reserved.
5// ORCID: 0009-0009-3560-0851
6// Contact: www.anulum.li | protoscience@anulum.li
7// SC-NeuroCore — Rust DNA circuit acceleration engine
8
9//! High-performance DNA strand displacement pipeline.
10//!
11//! Accelerates three critical hot paths from `sc_neurocore.bridges.dna_mapper`:
12//!
13//! 1. **Sequence design** — GC-balanced, homopolymer-free, orthogonal oligo generation
14//! 2. **Cross-hybridization** — O(n²) pairwise alignment scoring (rayon-parallelized)
15//! 3. **Kinetic simulation** — RK4 mass-action integrator with Arrhenius scaling
16
17pub(crate) mod bindings;
18
19use rand::RngExt;
20use rand::SeedableRng;
21use rand_xoshiro::Xoshiro256PlusPlus;
22use rayon::prelude::*;
23use std::collections::HashMap;
24
25// ── Constants ────────────────────────────────────────────────────────
26
27const ALPHABET: [u8; 4] = *b"ACGT";
28const GC_LOW: f64 = 0.40;
29const GC_HIGH: f64 = 0.60;
30const MAX_HOMOPOLYMER: usize = 3;
31#[allow(dead_code)] // reserved for strand-displacement toehold heuristic
32const DEFAULT_TOEHOLD: usize = 7;
33const R_GAS: f64 = 1.987e-3; // kcal/(mol·K)
34
35// ── Sequence Designer ────────────────────────────────────────────────
36
37/// Generate a GC-balanced, homopolymer-free DNA sequence.
38pub fn design_sequence(length: usize, seed: u64) -> Vec<u8> {
39    let mut rng = Xoshiro256PlusPlus::seed_from_u64(seed);
40    let mut seq = Vec::with_capacity(length);
41
42    for _ in 0..length * 10 {
43        if seq.len() >= length {
44            break;
45        }
46        let base = ALPHABET[rng.random_range(0..4)];
47
48        // Homopolymer check
49        if seq.len() >= MAX_HOMOPOLYMER {
50            let tail = &seq[seq.len() - MAX_HOMOPOLYMER..];
51            if tail.iter().all(|&b| b == base) {
52                continue;
53            }
54        }
55
56        seq.push(base);
57
58        // GC check at end
59        if seq.len() == length {
60            let gc = gc_content(&seq);
61            if !(GC_LOW..=GC_HIGH).contains(&gc) {
62                seq.clear();
63            }
64        }
65    }
66
67    // Fallback if constraints too tight
68    while seq.len() < length {
69        let base = ALPHABET[rng.random_range(0..4)];
70        seq.push(base);
71    }
72    seq.truncate(length);
73    seq
74}
75
76/// Generate multiple orthogonal sequences.
77pub fn design_orthogonal_set(count: usize, length: usize, seed: u64) -> Vec<Vec<u8>> {
78    let mut result = Vec::with_capacity(count);
79    for i in 0..count {
80        result.push(design_sequence(length, seed.wrapping_add(i as u64)));
81    }
82    result
83}
84
85#[inline]
86fn gc_content(seq: &[u8]) -> f64 {
87    if seq.is_empty() {
88        return 0.5;
89    }
90    let gc = seq.iter().filter(|&&b| b == b'G' || b == b'C').count();
91    gc as f64 / seq.len() as f64
92}
93
94/// Watson-Crick complement.
95pub fn complement(seq: &[u8]) -> Vec<u8> {
96    seq.iter()
97        .map(|&b| match b {
98            b'A' => b'T',
99            b'T' => b'A',
100            b'C' => b'G',
101            b'G' => b'C',
102            _ => b'N',
103        })
104        .collect()
105}
106
107/// Reverse complement.
108pub fn reverse_complement(seq: &[u8]) -> Vec<u8> {
109    let mut rc = complement(seq);
110    rc.reverse();
111    rc
112}
113
114// ── Cross-Hybridization Checker (rayon-parallelized) ─────────────────
115
116/// Pairwise alignment score between two sequences.
117fn alignment_score(a: &[u8], b: &[u8]) -> usize {
118    let rc_b = reverse_complement(b);
119    longest_common_substring(a, &rc_b)
120}
121
122fn longest_common_substring(a: &[u8], b: &[u8]) -> usize {
123    let n = a.len();
124    let m = b.len();
125    let mut max_len = 0usize;
126    let mut prev = vec![0usize; m + 1];
127    let mut curr = vec![0usize; m + 1];
128
129    for i in 1..=n {
130        for j in 1..=m {
131            if a[i - 1] == b[j - 1] {
132                curr[j] = prev[j - 1] + 1;
133                if curr[j] > max_len {
134                    max_len = curr[j];
135                }
136            } else {
137                curr[j] = 0;
138            }
139        }
140        std::mem::swap(&mut prev, &mut curr);
141        curr.iter_mut().for_each(|x| *x = 0);
142    }
143    max_len
144}
145
146/// Check all pairs for dangerous cross-hybridization.
147/// Returns vec of (i, j, score) for pairs above threshold.
148pub fn check_cross_hybridization(
149    sequences: &[Vec<u8>],
150    threshold: usize,
151) -> Vec<(usize, usize, usize)> {
152    let n = sequences.len();
153    if n < 2 {
154        return vec![];
155    }
156
157    // Build pair indices
158    let pairs: Vec<(usize, usize)> = (0..n)
159        .flat_map(|i| (i + 1..n).map(move |j| (i, j)))
160        .collect();
161
162    // Parallel pairwise scoring
163    pairs
164        .par_iter()
165        .filter_map(|&(i, j)| {
166            let score = alignment_score(&sequences[i], &sequences[j]);
167            if score >= threshold {
168                Some((i, j, score))
169            } else {
170                None
171            }
172        })
173        .collect()
174}
175
176// ── Kinetic Simulator (RK4 + Arrhenius) ──────────────────────────────
177
178/// Gate types for kinetic simulation.
179#[derive(Clone, Copy, Debug)]
180pub enum DnaGateType {
181    And,
182    Or,
183    Not,
184    Threshold,
185    Mux,
186    Amplifier,
187    Buffer,
188    Nand,
189    Xor,
190}
191
192/// A gate in the DNA circuit.
193#[derive(Clone, Debug)]
194pub struct DnaGateSpec {
195    pub gate_type: DnaGateType,
196    pub input_names: Vec<String>,
197    pub output_name: String,
198    pub threshold: f64,
199    pub leak_rate: f64,
200}
201
202/// Kinetic simulation configuration.
203pub struct KineticConfig {
204    pub k_hyb: f64,
205    pub k_disp: f64,
206    pub temperature_c: f64,
207    pub max_conc: f64,
208    pub use_rk4: bool,
209}
210
211impl Default for KineticConfig {
212    fn default() -> Self {
213        Self {
214            k_hyb: 3e5,
215            k_disp: 1.0,
216            temperature_c: 37.0,
217            max_conc: 200.0,
218            use_rk4: true,
219        }
220    }
221}
222
223/// Arrhenius temperature scaling.
224#[inline]
225fn arrhenius_scale(k_ref: f64, temperature_c: f64, ea_kcal: f64) -> f64 {
226    let t_ref = 310.15; // 37°C
227    let t_op = temperature_c + 273.15;
228    k_ref * (-(ea_kcal / R_GAS) * (1.0 / t_op - 1.0 / t_ref)).exp()
229}
230
231/// Compute effective rate constant for a gate.
232fn compute_k_eff(gate: &DnaGateSpec, inputs: &HashMap<String, f64>, config: &KineticConfig) -> f64 {
233    let k_hyb = arrhenius_scale(config.k_hyb, config.temperature_c, 15.0);
234    let k_disp = arrhenius_scale(config.k_disp, config.temperature_c, 15.0);
235
236    let k_eff = match gate.gate_type {
237        DnaGateType::And => {
238            let concs: Vec<f64> = gate
239                .input_names
240                .iter()
241                .map(|n| *inputs.get(n).unwrap_or(&0.0))
242                .collect();
243            let all_present = concs.iter().all(|&c| c > 0.0);
244            let min_c = concs.iter().cloned().fold(f64::MAX, f64::min);
245            k_hyb * min_c * 1e-9 * if all_present { 1.0 } else { 0.0 }
246        }
247        DnaGateType::Or => {
248            let concs: Vec<f64> = gate
249                .input_names
250                .iter()
251                .map(|n| *inputs.get(n).unwrap_or(&0.0))
252                .collect();
253            let max_c = concs.iter().cloned().fold(0.0f64, f64::max);
254            k_hyb * max_c * 1e-9
255        }
256        DnaGateType::Not => {
257            let inp = *inputs.get(&gate.input_names[0]).unwrap_or(&0.0);
258            k_disp * (1.0 - (inp / config.max_conc).min(1.0))
259        }
260        DnaGateType::Threshold => {
261            let inp = *inputs.get(&gate.input_names[0]).unwrap_or(&0.0);
262            let excess = (inp - gate.threshold * config.max_conc).max(0.0);
263            k_hyb * excess * 1e-9
264        }
265        DnaGateType::Mux => {
266            let sel = *inputs.get(&gate.input_names[0]).unwrap_or(&0.0);
267            let a = *inputs.get(&gate.input_names[1]).unwrap_or(&0.0);
268            let b = *inputs.get(&gate.input_names[2]).unwrap_or(&0.0);
269            let sel_frac = (sel / config.max_conc).min(1.0);
270            k_hyb * (sel_frac * a + (1.0 - sel_frac) * b) * 1e-9
271        }
272        DnaGateType::Amplifier => {
273            let inp = *inputs.get(&gate.input_names[0]).unwrap_or(&0.0);
274            k_hyb * inp * 1e-9 * 5.0
275        }
276        DnaGateType::Buffer => {
277            let inp = *inputs.get(&gate.input_names[0]).unwrap_or(&0.0);
278            k_disp * (inp / config.max_conc).min(1.0)
279        }
280        _ => 0.0,
281    };
282
283    k_eff + gate.leak_rate
284}
285
286/// Run kinetic simulation, returning time traces per output.
287pub fn simulate_kinetics(
288    gates: &[DnaGateSpec],
289    input_concentrations: &HashMap<String, f64>,
290    duration_s: f64,
291    dt: f64,
292    config: &KineticConfig,
293) -> HashMap<String, Vec<f64>> {
294    let n_steps = (duration_s / dt) as usize;
295    let max_conc = config.max_conc;
296
297    let mut result: HashMap<String, Vec<f64>> = HashMap::new();
298
299    // Time axis
300    let time: Vec<f64> = (0..n_steps).map(|t| t as f64 * dt).collect();
301    result.insert("time".to_string(), time);
302
303    for gate in gates {
304        let k_eff = compute_k_eff(gate, input_concentrations, config);
305        let mut conc = vec![0.0f64; n_steps];
306
307        if config.use_rk4 {
308            for t in 1..n_steps {
309                let c = conc[t - 1];
310                let k1 = k_eff * (max_conc - c) * dt;
311                let k2 = k_eff * (max_conc - (c + k1 / 2.0)) * dt;
312                let k3 = k_eff * (max_conc - (c + k2 / 2.0)) * dt;
313                let k4 = k_eff * (max_conc - (c + k3)) * dt;
314                conc[t] = (c + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0)
315                    .max(0.0)
316                    .min(max_conc);
317            }
318        } else {
319            for t in 1..n_steps {
320                let d = k_eff * (max_conc - conc[t - 1]) * dt;
321                conc[t] = (conc[t - 1] + d).max(0.0).min(max_conc);
322            }
323        }
324
325        result.insert(gate.output_name.clone(), conc);
326    }
327
328    result
329}
330
331// ── Hairpin Detection ────────────────────────────────────────────────
332
333/// Detect hairpins in a sequence. Returns vec of (stem_start, stem_len, loop_len).
334pub fn detect_hairpins(seq: &[u8], min_stem: usize, min_loop: usize) -> Vec<(usize, usize, usize)> {
335    let n = seq.len();
336    let mut hairpins = Vec::new();
337
338    if n < min_stem * 2 + min_loop {
339        return hairpins;
340    }
341
342    let wc = |a: u8, b: u8| -> bool {
343        matches!(
344            (a, b),
345            (b'A', b'T') | (b'T', b'A') | (b'C', b'G') | (b'G', b'C')
346        )
347    };
348
349    for i in 0..n.saturating_sub(min_stem * 2 + min_loop) {
350        for stem_len in min_stem..12.min((n - i) / 2) {
351            let loop_start = i + stem_len;
352            for loop_len in min_loop..10.min(n - loop_start - stem_len + 1) {
353                let j = loop_start + loop_len;
354                if j + stem_len > n {
355                    break;
356                }
357                let matches = (0..stem_len)
358                    .filter(|&k| wc(seq[i + k], seq[j + stem_len - 1 - k]))
359                    .count();
360                if matches >= stem_len {
361                    hairpins.push((i, stem_len, loop_len));
362                }
363            }
364        }
365    }
366    hairpins
367}
368
369// ── Tests ────────────────────────────────────────────────────────────
370
371#[cfg(test)]
372mod tests {
373    use super::*;
374
375    #[test]
376    fn test_design_sequence_length() {
377        let seq = design_sequence(30, 42);
378        assert_eq!(seq.len(), 30);
379    }
380
381    #[test]
382    fn test_design_sequence_alphabet() {
383        let seq = design_sequence(50, 42);
384        assert!(seq.iter().all(|&b| ALPHABET.contains(&b)));
385    }
386
387    #[test]
388    fn test_gc_content_balanced() {
389        let seq = design_sequence(40, 42);
390        let gc = gc_content(&seq);
391        assert!((0.3..=0.7).contains(&gc), "GC={gc}");
392    }
393
394    #[test]
395    fn test_complement() {
396        assert_eq!(complement(b"ACGT"), b"TGCA");
397    }
398
399    #[test]
400    fn test_reverse_complement() {
401        assert_eq!(reverse_complement(b"ACGT"), b"ACGT");
402    }
403
404    #[test]
405    fn test_cross_hybridization_self() {
406        let seqs = vec![b"ACGTACGTACGT".to_vec(), b"ACGTACGTACGT".to_vec()];
407        let flags = check_cross_hybridization(&seqs, 4);
408        assert!(!flags.is_empty());
409    }
410
411    #[test]
412    fn test_cross_hybridization_orthogonal() {
413        let seqs = design_orthogonal_set(3, 30, 42);
414        let flags = check_cross_hybridization(&seqs, 20);
415        assert!(flags.is_empty());
416    }
417
418    #[test]
419    fn test_simulate_kinetics_and_gate() {
420        let gates = vec![DnaGateSpec {
421            gate_type: DnaGateType::And,
422            input_names: vec!["A".to_string(), "B".to_string()],
423            output_name: "C".to_string(),
424            threshold: 0.5,
425            leak_rate: 1e-7,
426        }];
427        let mut inputs = HashMap::new();
428        inputs.insert("A".to_string(), 200.0);
429        inputs.insert("B".to_string(), 200.0);
430
431        let result = simulate_kinetics(&gates, &inputs, 1800.0, 1.0, &KineticConfig::default());
432        let c = result.get("C").unwrap();
433        assert!(c.last().unwrap() > &50.0);
434    }
435
436    #[test]
437    fn test_arrhenius_higher_temp_faster() {
438        let k37 = arrhenius_scale(3e5, 37.0, 15.0);
439        let k25 = arrhenius_scale(3e5, 25.0, 15.0);
440        assert!(k37 > k25);
441    }
442
443    #[test]
444    fn test_detect_hairpins_short() {
445        let hps = detect_hairpins(b"ACGT", 4, 3);
446        assert!(hps.is_empty());
447    }
448
449    #[test]
450    fn test_lcs() {
451        let score = longest_common_substring(b"ACGTACGT", b"ACGTACGT");
452        assert_eq!(score, 8);
453    }
454}
sc_neurocore_engine/dna.rs

sc_neurocore_engine/
dna.rs
