Nanoscale Self-Assembly: Programming Matter at the Molecular Scale

import numpy as np
from dataclasses import dataclass
from typing import List, Tuple

@dataclass
class DNAStrand:
    """Single DNA strand with sequence"""
    name: str
    sequence: str
    
    def complement(self):
        """Generate Watson-Crick complement"""
        comp_map = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G'}
        return ''.join(comp_map[base] for base in self.sequence[::-1])
    
    def binding_energy(self, other_seq: str, overlap: int) -> float:
        """Calculate hybridization energy (simplified)"""
        # AT bond: 2 kcal/mol, GC bond: 3 kcal/mol
        energy = 0.0
        for i in range(overlap):
            if i >= len(self.sequence) or i >= len(other_seq):
                break
            
            base1 = self.sequence[i]
            base2 = other_seq[i]
            
            if (base1 == 'A' and base2 == 'T') or (base1 == 'T' and base2 == 'A'):
                energy += -2.0
            elif (base1 == 'G' and base2 == 'C') or (base1 == 'C' and base2 == 'G'):
                energy += -3.0
        
        return energy

class DNAOrigami:
    """Design DNA origami nanostructures"""
    def __init__(self, scaffold_length=7249):  # M13mp18 scaffold
        self.scaffold_length = scaffold_length
        self.staples = []
        self.geometry = []
    
    def add_staple(self, start_pos: int, end_pos: int, crossover_pos: int = None):
        """Add staple strand that binds to scaffold"""
        length = abs(end_pos - start_pos)
        
        # Generate staple sequence (complement to scaffold region)
        # In practice, use actual M13mp18 scaffold sequence
        staple_seq = self._generate_staple_sequence(start_pos, end_pos)
        
        staple = {
            'start': start_pos,
            'end': end_pos,
            'sequence': staple_seq,
            'crossover': crossover_pos  # Connects to adjacent helix
        }
        
        self.staples.append(staple)
        
        # Validate staple design
        if length < 18:
            print(f"⚠️ Warning: Staple too short ({length} bp), may not bind stably")
        if length > 60:
            print(f"⚠️ Warning: Staple too long ({length} bp), may form secondary structure")
    
    def _generate_staple_sequence(self, start: int, end: int) -> str:
        """Generate staple complementary to scaffold region"""
        # Simplified: random sequence (real implementation uses actual scaffold)
        bases = ['A', 'T', 'G', 'C']
        length = abs(end - start)
        return ''.join(np.random.choice(bases) for _ in range(length))
    
    def design_rectangle(self, width_helices: int, length_bp: int):
        """Design simple rectangular origami"""
        print(f"Designing {width_helices}-helix × {length_bp}-bp rectangle\n")
        
        staple_length = 32  # Typical staple length
        
        # Place staples along each helix
        for helix_idx in range(width_helices):
            offset = helix_idx * self.scaffold_length // width_helices
            
            for pos in range(0, length_bp, staple_length):
                start = offset + pos
                end = offset + pos + staple_length
                
                # Crossover every 7 bp to adjacent helix (for rigidity)
                crossover = None
                if pos % 21 == 0 and helix_idx < width_helices - 1:
                    crossover = helix_idx + 1
                
                self.add_staple(start, end, crossover)
        
        print(f"Generated {len(self.staples)} staple strands")
    
    def export_sequences(self, filename: str):
        """Export staple sequences for DNA synthesis"""
        with open(filename, 'w') as f:
            f.write("Staple_Name,Sequence,Length\n")
            for i, staple in enumerate(self.staples):
                f.write(f"Staple_{i},{staple['sequence']},{len(staple['sequence'])}\n")
        
        print(f"Exported {len(self.staples)} sequences to {filename}")
        print(f"Total DNA synthesis cost: ~${len(self.staples) * 0.10:.2f} (at $0.10/staple)")

# Design DNA origami rectangle
origami = DNAOrigami()
origami.design_rectangle(width_helices=6, length_bp=100)
origami.export_sequences("staples.csv")

class MolecularReplicator:
    """
    ⚠️ WARNING: Self-replicating systems can grow exponentially
    This is a THEORETICAL model for educational purposes only
    """
    def __init__(self, replication_rate=1.0, error_rate=0.001):
        self.population = 1
        self.generation = 0
        self.replication_rate = replication_rate  # copies per timestep
        self.error_rate = error_rate  # mutation probability
        self.mutants = 0
        
        # Safety bounds
        self.max_population = 1e12  # Kill switch threshold
        self.max_generations = 100
    
    def replicate(self, timesteps=10):
        """Simulate exponential replication"""
        history = [self.population]
        
        for t in range(timesteps):
            # Exponential growth
            new_copies = int(self.population * self.replication_rate)
            
            # Mutations during replication
            new_mutants = int(new_copies * self.error_rate)
            self.mutants += new_mutants
            
            self.population += new_copies
            self.generation += 1
            
            history.append(self.population)
            
            # ⚠️ Safety checks
            if self.population > self.max_population:
                print(f"⚠️ CRITICAL: Population exceeded safe threshold at t={t}")
                print(f"   Population: {self.population:.2e}")
                print(f"   Mutants: {self.mutants}")
                print("   EMERGENCY SHUTDOWN TRIGGERED")
                break
            
            if self.generation > self.max_generations:
                print(f"⚠️ Max generations reached, halting replication")
                break
        
        return history
    
    def calculate_doubling_time(self):
        """Time for population to double"""
        return np.log(2) / np.log(1 + self.replication_rate)
    
    def estimate_resource_consumption(self, molecule_mass_g=1e-15):
        """Estimate mass consumed by replicators"""
        total_mass_g = self.population * molecule_mass_g
        
        # Convert to grams
        if total_mass_g < 1e-6:
            return f"{total_mass_g * 1e9:.2f} nanograms"
        elif total_mass_g < 1e-3:
            return f"{total_mass_g * 1e6:.2f} micrograms"
        elif total_mass_g < 1:
            return f"{total_mass_g * 1e3:.2f} milligrams"
        else:
            return f"{total_mass_g:.2f} grams"

# ⚠️ Demonstrate exponential growth danger
print("=== Self-Replicating Molecular System ===\n")
print("⚠️ WARNING: This simulation demonstrates why uncontrolled")
print("   self-replication is catastrophic\n")

replicator = MolecularReplicator(replication_rate=2.0, error_rate=0.01)

print(f"Initial population: {replicator.population}")
print(f"Replication rate: {replicator.replication_rate}x per generation")
print(f"Doubling time: {replicator.calculate_doubling_time():.2f} generations\n")

history = replicator.replicate(timesteps=30)

print(f"\nFinal population: {replicator.population:.2e}")
print(f"Generations: {replicator.generation}")
print(f"Mutant fraction: {replicator.mutants / replicator.population:.2%}")
print(f"Mass consumed: {replicator.estimate_resource_consumption()}")

# Grey goo scenario calculation
print("\n=== Grey Goo Calculation ===")
print("If replicators convert all carbon on Earth:")
earth_carbon_kg = 1.9e20  # kg of carbon
molecule_mass_g = 1e-15
max_replicators = earth_carbon_kg * 1000 / molecule_mass_g
generations_to_consume_earth = np.log2(max_replicators)
print(f"Replicators at saturation: {max_replicators:.2e}")
print(f"Generations to consume biosphere: ~{generations_to_consume_earth:.0f}")
print(f"At 1 generation/hour: ~{generations_to_consume_earth/24:.1f} days")
print("\n⚠️ This is why containment is CRITICAL")

class DNATile:
    """DNA tile with sticky ends for programmed assembly"""
    def __init__(self, name: str, north: str, east: str, south: str, west: str):
        self.name = name
        # Sticky ends (4-letter codes)
        self.edges = {
            'N': north,
            'E': east,
            'S': south,
            'W': west
        }
    
    def can_bind(self, other: 'DNATile', direction: str) -> bool:
        """Check if tiles can bind in given direction"""
        opposite = {'N': 'S', 'S': 'N', 'E': 'W', 'W': 'E'}
        
        my_edge = self.edges[direction]
        other_edge = other.edges[opposite[direction]]
        
        # Complementary sticky ends bind
        return my_edge == complement(other_edge)

def complement(seq: str) -> str:
    """Watson-Crick complement"""
    comp = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G'}
    return ''.join(comp[b] for b in seq[::-1])

class TileAssembly:
    """Simulate 2D DNA tile self-assembly"""
    def __init__(self, tile_types: List[DNATile], seed_tile: DNATile):
        self.tile_types = tile_types
        self.grid = {(0, 0): seed_tile}  # Start from seed
        self.assembly_steps = 0
    
    def grow(self, max_steps=1000):
        """Grow crystal by adding compatible tiles"""
        for step in range(max_steps):
            # Find all growth sites (empty neighbors of existing tiles)
            growth_sites = set()
            for (x, y) in self.grid.keys():
                for dx, dy, direction in [(0, 1, 'N'), (1, 0, 'E'), (0, -1, 'S'), (-1, 0, 'W')]:
                    neighbor = (x + dx, y + dy)
                    if neighbor not in self.grid:
                        growth_sites.add((neighbor, direction))
            
            if not growth_sites:
                break  # No more growth sites
            
            # Try to add a tile at each site
            for (x, y), direction in growth_sites:
                # Find compatible tile
                for tile in self.tile_types:
                    # Check binding to existing neighbors
                    compatible = True
                    for dx, dy, dir_check in [(0, 1, 'N'), (1, 0, 'E'), (0, -1, 'S'), (-1, 0, 'W')]:
                        neighbor_pos = (x + dx, y + dy)
                        if neighbor_pos in self.grid:
                            neighbor_tile = self.grid[neighbor_pos]
                            if not tile.can_bind(neighbor_tile, dir_check):
                                compatible = False
                                break
                    
                    if compatible:
                        self.grid[(x, y)] = tile
                        self.assembly_steps += 1
                        break
        
        return self.grid

# Example: Binary counter using DNA tiles
tile_0 = DNATile("Zero", north="ATCG", east="GCTA", south="CGAT", west="TACG")
tile_1 = DNATile("One", north="TAGC", east="CGTA", south="ATGC", west="GCTA")

assembly = TileAssembly([tile_0, tile_1], seed_tile=tile_0)
final_structure = assembly.grow(max_steps=100)

print(f"Assembly complete: {len(final_structure)} tiles, {assembly.assembly_steps} steps")