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

sc_neurocore_engine/
dna.rs
