fiona - More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc

I'm attempting to do a spatial join much like the example here: Is there a python option to "join attributes by location"?. However, that approach seems really inefficient / slow. Even running this with a modest 250 points takes almost 2 minutes and it fails entirely on shapefiles with > 1,000 points. Is there a better approach? I'd like to do this entirely in Python without using ArcGIS, QGIS, etc.

I'd also be interested to know if it's possible to SUM attributes (i.e. population) of all the points that fall within a polygon and join that quantity to the polygon shapefile.

Here is the code I'm trying to convert. I get an error on line 9:

poly['properties']['score'] += point['properties']['score']

which says:

TypeError: unsupported operand type(s) for +=: 'NoneType' and 'float'.

If I replace the "+=" with "=" it runs fine but that doesn't sum the fields. I've also tried making these as integers but that fails as well.

with fiona.open(poly_shp, 'r') as n: 
  with fiona.open(point_shp,'r') as s:

    outSchema = {'geometry': 'Polygon','properties':{'region':'str','score':'float'}}
    with fiona.open (out_shp, 'w', 'ESRI Shapefile', outSchema, crs) as output:
        for point in s:
            for poly in n:
                if shape(point['geometry']).within(shape(poly['geometry'])):  
                    poly['properties']['score']) += point['properties']['score'])
                    output.write({
                        'properties':{
                            'region':poly['properties']['NAME'],
                            'score':poly['properties']['score']},

                        'geometry':poly['geometry']})