Simulating cell population dynamics with Haskell

telo.hs

data Cell = Cell {teloLength :: Double, birthDate :: Double, lifeSpan :: Double, seedCell :: Integer}
instance Eq Cell where (Cell t1 _ _ _) == (Cell t2 _ _ _) = t1 == t2
instance Ord Cell where (Cell t1 _ _ _) `compare` (Cell t2 _ _ _) = t1 `compare` t2
instance Show Cell where show (Cell t b _ _) = show t ++ " " ++ show b

data CellLine x = Senescent | Mother x (CellLine x) (CellLine x) deriving (Show, Read, Eq) 

-- Functions modelling cell biology
shortenTelo :: Cell -> Cell
shortenTelo (Cell telo x y oldseed) = Cell (telo - delta) x y newseed
  where 
  newseed = nextRand oldseed
  u1 = makeUnif oldseed
  u2 = makeUnif newseed
  delta = boxMuller 1 0.1 (u1,u2)
  
divideCell :: Cell -> CellLine Cell
divideCell (Cell telo bday lifespan seed)
	| telo <= 0 = Senescent
	| telo > 0  = Mother (Cell telo bday lifespan seed) (divideCell (shortenTelo (Cell telo (bday+lifespan) delta seed2))) (divideCell (shortenTelo (Cell telo (bday+lifespan) delta seed3)))
	where 
	seed1 = nextRand seed
	seed2 = nextRand seed1
	seed3 = nextRand seed2
	delta = boxMuller 1 0.2 (makeUnif seed, makeUnif seed1) 

-- Functions simulating measurements
maxDivisions :: CellLine x -> Integer
maxDivisions Senescent = 1
maxDivisions (Mother x daughterLeft daughterRight) = 1 + max (maxDivisions daughterLeft) (maxDivisions daughterRight)

countCells :: CellLine Cell -> Double -> Integer
countCells Senescent x = 0
countCells (Mother m left right) x = if ((birthDate$m) < x) && (x <= ((birthDate$m) + (lifeSpan$m)))
	then 1 + countCells (left) x + countCells (right) x
	else countCells (left) x + countCells (right) x

countLeaves :: CellLine Cell -> Integer
countLeaves Senescent = 1
countLeaves (Mother x left right)
  | x <= (Cell 0 0 0 0) = 1 + countLeaves(left) + countLeaves(right)
  | x > (Cell 0 0 0 0) = countLeaves(left) + countLeaves(right)

-- Functions & data for generating sequences of pseudo-random numbers
seed=17489
multiplier=25173
increment=13849
modulus=65536

nextRand :: Integer -> Integer
nextRand n = (multiplier*n + increment) `mod` modulus

makeUnif :: Integer -> Double
makeUnif seed = (fromIntegral seed :: Double)/(fromIntegral(modulus-1) :: Double)

boxMuller :: Double -> Double -> (Double, Double) -> Double
boxMuller mu sigma (r1,r2) =  mu + sigma * sqrt (-2 * log r1) * cos (2 * pi * r2)

-- Simulate
clin = divideCell (Cell 4 0 1 123)
mdivs = maxDivisions clin
numcells = countCells clin 2
numleaves = countLeaves clin