Updating raster value to reflect total number of points within cell using Python?

Using any of the tools you mentioned you can preprocess the points so that one feature (with a count) occurs in each raster cell. Converting that to a raster using the standard built-in method completes the task.

The advantage of this method is that it provides complete control over the conversion process, so you know exactly what is happening. Moreover, it easily generalizes: with only trivial modifications (at the summarize step) you can produce rasters whose values are any statistical summary of the points falling within each cell (of which the count is a simple example).

The preprocessing is relatively simple: you need to identify the cell in which any point (x,y) falls. The raster's extent is given by an origin (Ox, Oy) and corresponding cell sizes (h, k) (usually h=k). The (column, row) coordinates of (x,y)'s cell are

(i,j) = floor( (x-Ox)/h, (y-Oy)/k )

That requires two quick calculations to create new shapefile attributes i and j. Summarize the attribute table on the concatenation of i and j (which will identify the cells), retaining the total count (or any other summary) for each unique (i,j) pair, along with the values of i and j themselves. Convert the (i,j) attributes in the summary table back to coordinates of cell centers as

(x,y) = ( Ox + h*(i+1/2), Oy + k*(j+1/2) )

This produces a table with three attributes: x, y, and the summary value for the points corresponding to (x,y)'s cell. Proceed in any way you prefer to convert this into a raster.

I suggest using the Point to Raster tool. With default settings, it chooses one point per cell to generate raster values. Using the cell_assignment option of COUNT, however, the raster stores how many points occurred in each cell.

This produces new rasters rather than working with existing ones, but the process of combining rasters for analysis is simpler than combining raster/vector data.

Example arcpy code (see help link above for full context) would be:

# Set local variables
inFeatures = "ca_ozone_pts.shp"
valField = "ELEVATION"
outRaster = "c:/output/ca_elev02"
assignmentType = "COUNT"
priorityField = ""
cellSize = 2000

# Execute PointToRaster
arcpy.PointToRaster_conversion(inFeatures, valField, outRaster, 
                               assignmentType, priorityField, cellSize)